diff --git a/benchmarks/linear_programming/cuopt/initial_problem_check.hpp b/benchmarks/linear_programming/cuopt/initial_problem_check.hpp index fb89d24f0..73ff4f1fa 100644 --- a/benchmarks/linear_programming/cuopt/initial_problem_check.hpp +++ b/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& op_problem, - const std::vector& 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& op_problem, + const std::vector& primal_vars, + double abs_tol, + double rel_tol, + double int_tol = 1e-5) { const std::vector& values = op_problem.get_constraint_matrix_values(); const std::vector& indices = op_problem.get_constraint_matrix_indices(); @@ -52,6 +53,7 @@ void test_constraint_sanity(const cuopt::mps_parser::mps_data_model_t& constraint_upper_bounds = op_problem.get_constraint_upper_bounds(); const std::vector& variable_lower_bounds = op_problem.get_variable_lower_bounds(); const std::vector& variable_upper_bounds = op_problem.get_variable_upper_bounds(); + const std::vector& variable_types = op_problem.get_variable_types(); std::vector residual(constraint_lower_bounds.size(), 0.0); std::vector viol(constraint_lower_bounds.size(), 0.0); @@ -64,6 +66,7 @@ void test_constraint_sanity(const cuopt::mps_parser::mps_data_model_t{}; + 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(constraint_lower_bounds[i], @@ -71,25 +74,33 @@ void test_constraint_sanity(const cuopt::mps_parser::mps_data_model_t 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; } \ No newline at end of file diff --git a/benchmarks/linear_programming/cuopt/initial_solution_reader.hpp b/benchmarks/linear_programming/cuopt/initial_solution_reader.hpp index a8f9ee564..85d2b80ef 100644 --- a/benchmarks/linear_programming/cuopt/initial_solution_reader.hpp +++ b/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 #include #include @@ -24,7 +25,7 @@ class solution_reader_t { public: std::unordered_map 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); diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index c24b8b39c..f3124136f 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/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(); } -// reads a solution from an input file. The input file needs to be csv formatted -// var_name,val -std::vector read_solution_from_file(const std::string file_path, - const std::vector& var_names) +void read_single_solution_from_path(const std::string& path, + const std::vector& var_names, + std::vector>& 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 assignment; for (auto name : var_names) { @@ -112,20 +111,42 @@ std::vector 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> read_solution_from_dir(const std::string file_path, + const std::string& mps_file_name, + const std::vector& var_names) +{ + std::vector> 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 input_file_dir, + std::optional 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; diff --git a/cpp/cuopt_cli.cpp b/cpp/cuopt_cli.cpp index 594fb745d..7394dc7e7 100644 --- a/cpp/cuopt_cli.cpp +++ b/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 { diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 24257f920..fbb75d80b 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/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::max(); +}; + // Forward declare solver_settings_t for friend class template class solver_settings_t; @@ -44,26 +50,19 @@ 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 - */ - rmm::device_uvector& get_initial_solution() const; - /** * @brief Get the callback for the user solution * @@ -71,8 +70,6 @@ class mip_solver_settings_t { */ const std::vector 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> initial_solution_; + /** Initial primal solutions */ + std::vector>> 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 mip_callbacks_; diff --git a/cpp/include/cuopt/linear_programming/solver_settings.hpp b/cpp/include/cuopt/linear_programming/solver_settings.hpp index 5512b9de1..f5b223183 100644 --- a/cpp/include/cuopt/linear_programming/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/solver_settings.hpp @@ -88,10 +88,11 @@ class solver_settings_t { const rmm::device_uvector& 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& get_initial_mip_solution() const; const pdlp_warm_start_data_view_t& get_pdlp_warm_start_data_view() const noexcept; const std::vector get_mip_callbacks() const; diff --git a/cpp/src/linear_programming/utilities/problem_checking.cu b/cpp/src/linear_programming/utilities/problem_checking.cu index b7104fd8f..57f10e55b 100644 --- a/cpp/src/linear_programming/utilities/problem_checking.cu +++ b/cpp/src/linear_programming/utilities/problem_checking.cu @@ -113,8 +113,8 @@ void problem_checking_t::check_initial_solution_representation( const optimization_problem_t& op_problem, const mip_solver_settings_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); } } diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index ac60c168c..0e148610c 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -371,9 +371,11 @@ const rmm::device_uvector& solver_settings_t::get_initial_pdlp_du } template -void solver_settings_t::set_initial_mip_solution(const f_t* solution, i_t size) +void solver_settings_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 @@ -382,12 +384,6 @@ void solver_settings_t::set_mip_callback(internals::base_solution_call mip_settings.set_mip_callback(callback); } -template -const rmm::device_uvector& solver_settings_t::get_initial_mip_solution() const -{ - return mip_settings.get_initial_solution(); -} - template const std::vector solver_settings_t::get_mip_callbacks() const diff --git a/cpp/src/mip/diversity/diversity_config.hpp b/cpp/src/mip/diversity/diversity_config.hpp new file mode 100644 index 000000000..30e7d25e4 --- /dev/null +++ b/cpp/src/mip/diversity/diversity_config.hpp @@ -0,0 +1,45 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +namespace cuopt::linear_programming::detail { + +struct diversity_config_t { + static constexpr double time_ratio_on_init_lp = 0.1; + static constexpr double max_time_on_lp = 30; + static constexpr double time_ratio_of_probing_cache = 0.10; + static constexpr double max_time_on_probing = 60; + static constexpr size_t max_iterations_without_improvement = 15; + static constexpr int max_var_diff = 256; + static constexpr size_t max_solutions = 32; + static constexpr double initial_infeasibility_weight = 1000.; + static constexpr double default_time_limit = 10.; + static constexpr int initial_island_size = 3; + static constexpr int maximum_island_size = 8; + static constexpr bool use_avg_diversity = false; + static constexpr double generation_time_limit_ratio = 0.6; + static constexpr double max_island_gen_time = 600; + static constexpr size_t n_sol_for_skip_init_gen = 3; + static constexpr double max_fast_sol_time = 10; + static constexpr double lp_run_time_if_feasible = 15.; + static constexpr double lp_run_time_if_infeasible = 1; + static constexpr double close_to_parents_ratio = 0.1; + static constexpr bool halve_population = false; +}; + +} // namespace cuopt::linear_programming::detail \ No newline at end of file diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu index f04a2a814..874230a29 100644 --- a/cpp/src/mip/diversity/diversity_manager.cu +++ b/cpp/src/mip/diversity/diversity_manager.cu @@ -16,25 +16,35 @@ */ #include -#include -#include "diversity_manager.cuh" - #include #include +#include +#include "diversity_config.hpp" +#include "diversity_manager.cuh" #include #include "cuda_profiler_api.h" +constexpr bool from_dir = false; +constexpr bool fj_only_run = false; + namespace cuopt::linear_programming::detail { -constexpr int max_var_diff = 256; -constexpr size_t max_solutions = 32; -constexpr double initial_infeasibility_weight = 1000.; -constexpr double default_time_limit = 10.; -constexpr int initial_island_size = 3; -constexpr int maximum_island_size = 8; -constexpr bool use_avg_diversity = false; +constexpr int max_var_diff = diversity_config_t::max_var_diff; +constexpr size_t max_solutions = diversity_config_t::max_solutions; +constexpr double initial_infeasibility_weight = diversity_config_t::initial_infeasibility_weight; +constexpr double default_time_limit = diversity_config_t::default_time_limit; +constexpr int initial_island_size = diversity_config_t::initial_island_size; +constexpr int maximum_island_size = diversity_config_t::maximum_island_size; +constexpr bool use_avg_diversity = diversity_config_t::use_avg_diversity; + +size_t fp_recombiner_config_t::max_n_of_vars_from_other = + fp_recombiner_config_t::initial_n_of_vars_from_other; +size_t ls_recombiner_config_t::max_n_of_vars_from_other = + ls_recombiner_config_t::initial_n_of_vars_from_other; +size_t bp_recombiner_config_t::max_n_of_vars_from_other = + bp_recombiner_config_t::initial_n_of_vars_from_other; template diversity_manager_t::diversity_manager_t(mip_solver_context_t& context_) @@ -60,8 +70,37 @@ diversity_manager_t::diversity_manager_t(mip_solver_context_thandle_ptr), rng(cuopt::seed_generator::get_seed()), - stats(context.stats) + stats(context.stats), + mab_arm_stats_(recombiner_enum_t::SIZE), + mab_rng_(cuopt::seed_generator::get_seed()) { + // Read configuration ID from environment variable + int max_config = -1; + // Read max configuration value from environment variable + const char* env_max_config = std::getenv("CUOPT_MAX_CONFIG"); + if (env_max_config != nullptr) { + try { + max_config = std::stoi(env_max_config); + CUOPT_LOG_DEBUG("Using maximum configuration value from environment: %d", max_config); + } catch (const std::exception& e) { + CUOPT_LOG_WARN("Failed to parse CUOPT_MAX_CONFIG environment variable: %s", e.what()); + } + } + if (max_config > 1) { + int config_id = -1; // Default value + const char* env_config_id = std::getenv("CUOPT_CONFIG_ID"); + if (env_config_id != nullptr) { + try { + config_id = std::stoi(env_config_id); + CUOPT_LOG_DEBUG("Using configuration ID from environment: %d", config_id); + } catch (const std::exception& e) { + CUOPT_LOG_WARN("Failed to parse CUOPT_CONFIG_ID environment variable: %s", e.what()); + } + } + run_only_ls_recombiner = config_id == 0; + run_only_bp_recombiner = config_id == 1; + run_only_fp_recombiner = config_id == 2; + } } // There should be at least 3 solutions in the population @@ -94,7 +133,7 @@ std::vector> diversity_manager_t::generate_more_s timer_t total_time_to_generate = timer_t(timer.remaining_time() / 5.); f_t time_limit = std::min(60., total_time_to_generate.remaining_time()); f_t ls_limit = std::min(5., timer.remaining_time() / 20.); - const i_t n_sols_to_generate = 2; + const i_t n_sols_to_generate = 3; for (i_t i = 0; i < n_sols_to_generate; ++i) { CUOPT_LOG_DEBUG("Trying to generate more solutions"); time_limit = std::min(time_limit, timer.remaining_time()); @@ -143,13 +182,12 @@ void diversity_manager_t::average_fj_weights(i_t i) } template -void diversity_manager_t::add_user_given_solution( +void diversity_manager_t::add_user_given_solutions( std::vector>& initial_sol_vector) { - if (context.settings.has_initial_solution()) { + for (const auto& init_sol : context.settings.initial_solutions) { solution_t sol(*problem_ptr); - auto& init_sol = context.settings.get_initial_solution(); - rmm::device_uvector init_sol_assignment(init_sol, sol.handle_ptr->get_stream()); + rmm::device_uvector init_sol_assignment(*init_sol, sol.handle_ptr->get_stream()); if (problem_ptr->pre_process_assignment(init_sol_assignment)) { raft::copy(sol.assignment.data(), init_sol_assignment.data(), @@ -169,7 +207,7 @@ void diversity_manager_t::add_user_given_solution( Assignment size %lu \ initial solution size %lu", sol.assignment.size(), - init_sol.size()); + init_sol->size()); } } } @@ -179,16 +217,19 @@ void diversity_manager_t::add_user_given_solution( template void diversity_manager_t::generate_initial_solutions() { - add_user_given_solution(initial_sol_vector); + add_user_given_solutions(initial_sol_vector); + bool skip_initial_island_generation = + initial_sol_vector.size() > diversity_config_t::n_sol_for_skip_init_gen || from_dir; // allocate maximum of 40% of the time to the initial island generation // aim to generate at least 5 feasible solutions thus spending 8% of the time to generate a // solution if we can generate faster generate up to 10 sols - const f_t generation_time_limit = 0.6 * timer.get_time_limit(); - const f_t max_island_gen_time = 600; - f_t total_island_gen_time = std::min(generation_time_limit, max_island_gen_time); + const f_t generation_time_limit = + diversity_config_t::generation_time_limit_ratio * timer.get_time_limit(); + const f_t max_island_gen_time = diversity_config_t::max_island_gen_time; + f_t total_island_gen_time = std::min(generation_time_limit, max_island_gen_time); timer_t gen_timer(total_island_gen_time); f_t sol_time_limit = gen_timer.remaining_time(); - for (i_t i = 0; i < maximum_island_size; ++i) { + for (i_t i = 0; i < maximum_island_size && !skip_initial_island_generation; ++i) { if (check_b_b_preemption()) { return; } if (i + population.get_external_solution_size() >= 5) { break; } CUOPT_LOG_DEBUG("Generating sol %d", i); @@ -232,7 +273,7 @@ void diversity_manager_t::generate_initial_solutions() population.var_threshold); population.print(); auto new_sol_vector = population.get_external_solutions(); - recombine_and_ls_with_all(new_sol_vector); + if (!fj_only_run) { recombine_and_ls_with_all(new_sol_vector); } } template @@ -267,7 +308,8 @@ void diversity_manager_t::generate_quick_feasible_solution() { solution_t solution(*problem_ptr); // min 1 second, max 10 seconds - const f_t generate_fast_solution_time = std::min(10., std::max(1., timer.remaining_time() / 20.)); + const f_t generate_fast_solution_time = + std::min(diversity_config_t::max_fast_sol_time, std::max(1., timer.remaining_time() / 20.)); timer_t sol_timer(generate_fast_solution_time); // do very short LP run to get somewhere close to the optimal point ls.generate_fast_solution(solution, sol_timer); @@ -303,17 +345,31 @@ bool diversity_manager_t::check_b_b_preemption() // returns the best feasible solution template -solution_t diversity_manager_t::run_solver() +void diversity_manager_t::run_fj_alone(solution_t& solution) { - const f_t time_limit = timer.remaining_time(); - constexpr f_t time_ratio_on_init_lp = 0.1; - constexpr f_t max_time_on_lp = 30; - const f_t lp_time_limit = std::min(max_time_on_lp, time_limit * time_ratio_on_init_lp); + CUOPT_LOG_INFO("Running FJ alone!"); + solution.round_nearest(); + ls.fj.settings.mode = fj_mode_t::EXIT_NON_IMPROVING; + ls.fj.settings.n_of_minimums_for_exit = 20000 * 1000; + ls.fj.settings.update_weights = true; + ls.fj.settings.feasibility_run = false; + ls.fj.settings.termination = fj_termination_flags_t::FJ_TERMINATION_TIME_LIMIT; + ls.fj.settings.time_limit = timer.remaining_time(); + ls.fj.solve(solution); + CUOPT_LOG_INFO("FJ alone finished!"); +} +// returns the best feasible solution +template +solution_t diversity_manager_t::run_solver() +{ + population.timer = timer; + const f_t time_limit = timer.remaining_time(); + const f_t lp_time_limit = std::min(diversity_config_t::max_time_on_lp, + time_limit * diversity_config_t::time_ratio_on_init_lp); // to automatically compute the solving time on scope exit auto timer_raii_guard = cuopt::scope_guard([&]() { stats.total_solve_time = timer.elapsed_time(); }); - // after every change to the problem, we should resize all the relevant vars // we need to encapsulate that to prevent repetitions lp_optimal_solution.resize(problem_ptr->n_variables, problem_ptr->handle_ptr->get_stream()); @@ -322,6 +378,8 @@ solution_t diversity_manager_t::run_solver() ls.lb_constraint_prop.bounds_update.setup(ls.lb_constraint_prop.temp_problem); ls.constraint_prop.bounds_update.resize(*problem_ptr); problem_ptr->check_problem_representation(true); + // have the structure ready for reusing later + problem_ptr->compute_integer_fixed_problem(); // test problem is not ii cuopt_func_call( ls.constraint_prop.bounds_update.calculate_activity_on_problem_bounds(*problem_ptr)); @@ -331,53 +389,64 @@ solution_t diversity_manager_t::run_solver() population.initialize_population(); if (check_b_b_preemption()) { return population.best_feasible(); } // before probing cache or LP, run FJ to generate initial primal feasible solution - generate_quick_feasible_solution(); - constexpr f_t time_ratio_of_probing_cache = 0.10; - constexpr f_t max_time_on_probing = 60; + if (!from_dir) { generate_quick_feasible_solution(); } + constexpr f_t time_ratio_of_probing_cache = diversity_config_t::time_ratio_of_probing_cache; + constexpr f_t max_time_on_probing = diversity_config_t::max_time_on_probing; f_t time_for_probing_cache = std::min(max_time_on_probing, time_limit * time_ratio_of_probing_cache); timer_t probing_timer{time_for_probing_cache}; if (check_b_b_preemption()) { return population.best_feasible(); } - compute_probing_cache(ls.constraint_prop.bounds_update, *problem_ptr, probing_timer); + if (!fj_only_run) { + compute_probing_cache(ls.constraint_prop.bounds_update, *problem_ptr, probing_timer); + } // careful, assign the correct probing cache ls.lb_constraint_prop.bounds_update.probing_cache.probing_cache = ls.constraint_prop.bounds_update.probing_cache.probing_cache; if (check_b_b_preemption()) { return population.best_feasible(); } - lp_state_t& lp_state = context.lp_state; + lp_state_t& lp_state = problem_ptr->lp_state; // resize because some constructor might be called before the presolve lp_state.resize(*problem_ptr, problem_ptr->handle_ptr->get_stream()); - auto lp_result = get_relaxed_lp_solution(*problem_ptr, - lp_optimal_solution, - lp_state, - context.settings.tolerances.absolute_tolerance, - lp_time_limit, - false, - false); - population.allocate_solutions(); - ls.lp_optimal_exists = true; - if (lp_result.get_termination_status() == pdlp_termination_status_t::Optimal) { - // get lp user objective and pass it to set_new_user_bound - set_new_user_bound(lp_result.get_objective_value()); - } else if (lp_result.get_termination_status() == pdlp_termination_status_t::PrimalInfeasible) { - // PDLP's infeasibility detection isn't an exact method and might be subject to false positives. - // Issue a warning, and continue solving. - CUOPT_LOG_WARN("PDLP detected primal infeasibility, problem might be infeasible!"); - ls.lp_optimal_exists = false; - } else if (lp_result.get_termination_status() == pdlp_termination_status_t::DualInfeasible) { - CUOPT_LOG_WARN("PDLP detected dual infeasibility, problem might be unbounded!"); - ls.lp_optimal_exists = false; - } else if (lp_result.get_termination_status() == pdlp_termination_status_t::TimeLimit) { - CUOPT_LOG_DEBUG( - "Initial LP run exceeded time limit, continuing solver with partial LP result!"); - // note to developer, in debug mode the LP run might be too slow and it might cause PDLP not to - // bring variables within the bounds + relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = lp_time_limit; + lp_settings.tolerance = context.settings.tolerances.absolute_tolerance; + lp_settings.return_first_feasible = false; + lp_settings.save_state = true; + if (!fj_only_run) { + auto lp_result = + get_relaxed_lp_solution(*problem_ptr, lp_optimal_solution, lp_state, lp_settings); + ls.lp_optimal_exists = true; + if (lp_result.get_termination_status() == pdlp_termination_status_t::Optimal) { + set_new_user_bound(lp_result.get_objective_value()); + } else if (lp_result.get_termination_status() == pdlp_termination_status_t::PrimalInfeasible) { + CUOPT_LOG_ERROR("Problem is primal infeasible, continuing anyway!"); + ls.lp_optimal_exists = false; + } else if (lp_result.get_termination_status() == pdlp_termination_status_t::DualInfeasible) { + CUOPT_LOG_ERROR("PDLP detected dual infeasibility, continuing anyway!"); + ls.lp_optimal_exists = false; + } else if (lp_result.get_termination_status() == pdlp_termination_status_t::TimeLimit) { + CUOPT_LOG_DEBUG( + "Initial LP run exceeded time limit, continuing solver with partial LP result!"); + // note to developer, in debug mode the LP run might be too slow and it might cause PDLP not + // to bring variables within the bounds + } + // in case the pdlp returned var boudns that are out of bounds + clamp_within_var_bounds(lp_optimal_solution, problem_ptr, problem_ptr->handle_ptr); } - // in case the pdlp returned var boudns that are out of bounds - clamp_within_var_bounds(lp_optimal_solution, problem_ptr, problem_ptr->handle_ptr); + population.allocate_solutions(); if (check_b_b_preemption()) { return population.best_feasible(); } // generate a population with 5 solutions(FP+FJ) generate_initial_solutions(); + if (context.settings.benchmark_info_ptr != nullptr) { + context.settings.benchmark_info_ptr->objective_of_initial_population = + population.best_feasible().get_user_objective(); + } + + if (fj_only_run) { + run_fj_alone(population.best_feasible()); + return population.best_feasible(); + } + if (timer.check_time_limit()) { return population.best_feasible(); } main_loop(); return population.best_feasible(); @@ -387,8 +456,9 @@ template void diversity_manager_t::diversity_step() { // TODO when the solver is faster, increase this number - constexpr i_t max_iterations_without_improvement = 15; - bool improved = true; + constexpr i_t max_iterations_without_improvement = + diversity_config_t::max_iterations_without_improvement; + bool improved = true; while (improved) { int k = max_iterations_without_improvement; improved = false; @@ -396,7 +466,7 @@ void diversity_manager_t::diversity_step() if (check_b_b_preemption()) { return; } auto new_sol_vector = population.get_external_solutions(); recombine_and_ls_with_all(new_sol_vector); - + population.adjust_weights_according_to_best_feasible(); cuopt_assert(population.test_invariant(), ""); if (population.current_size() < 2) { CUOPT_LOG_DEBUG("Population degenerated in diversity step"); @@ -443,7 +513,6 @@ void diversity_manager_t::recombine_and_ls_with_all(solution_t(solution))); } template @@ -453,11 +522,16 @@ void diversity_manager_t::recombine_and_ls_with_all( raft::common::nvtx::range fun_scope("recombine_and_ls_with_all"); if (solutions.size() > 0) { CUOPT_LOG_INFO("Running recombiners on B&B solutions with size %lu", solutions.size()); + // add all solutions because time limit might have been consumed and we might have exited before for (auto& sol : solutions) { - population.add_solution(std::move(solution_t(sol))); cuopt_func_call(sol.test_feasibility(true)); + population.add_solution(std::move(solution_t(sol))); + } + for (auto& sol : solutions) { + if (timer.check_time_limit()) { return; } solution_t ls_solution(sol); ls.run_local_search(ls_solution, population.weights, timer); + if (timer.check_time_limit()) { return; } // TODO try if running LP with integers fixed makes it feasible if (ls_solution.get_feasible()) { CUOPT_LOG_DEBUG("External LS searched solution feasible, running recombiners!"); @@ -473,7 +547,7 @@ void diversity_manager_t::recombine_and_ls_with_all( template void diversity_manager_t::main_loop() { - population.diversity_start_time = timer.elapsed_time(); + population.start_threshold_adjustment(); recombine_stats.reset(); while (true) { if (check_b_b_preemption()) { break; } @@ -490,15 +564,26 @@ void diversity_manager_t::main_loop() if (timer.check_time_limit()) { break; } diversity_step(); if (timer.check_time_limit()) { break; } - population.halve_the_population(); - auto new_solutions = generate_more_solutions(); - auto current_population = population.population_to_vector(); - population.clear(); - current_population.insert(current_population.end(), - std::make_move_iterator(new_solutions.begin()), - std::make_move_iterator(new_solutions.end())); - population.find_diversity(current_population, use_avg_diversity); - population.add_solutions_from_vec(std::move(current_population)); + + if (diversity_config_t::halve_population) { + population.adjust_threshold(timer); + i_t prev_threshold = population.var_threshold; + population.halve_the_population(); + auto new_solutions = generate_more_solutions(); + auto current_population = population.population_to_vector(); + population.clear(); + current_population.insert(current_population.end(), + std::make_move_iterator(new_solutions.begin()), + std::make_move_iterator(new_solutions.end())); + population.find_diversity(current_population, use_avg_diversity); + // if the threshold is lower than the threshold we progress with time + // set it to the higher threshold + // population.var_threshold = max(population.var_threshold, prev_threshold); + population.add_solutions_from_vec(std::move(current_population)); + } else { + // increase the threshold/decrease the diversity + population.adjust_threshold(timer); + } // population.add_solutions_from_vec(std::move(new_solutions)); // idea to try, we can average the weights of the new solutions population.update_weights(); @@ -510,6 +595,28 @@ void diversity_manager_t::main_loop() population.print(); } +template +void diversity_manager_t::check_better_than_both(solution_t& offspring, + solution_t& sol1, + solution_t& sol2) +{ + bool better_than_both = false; + if (sol1.get_feasible() && sol2.get_feasible()) { + better_than_both = offspring.get_objective() < + (std::min(sol1.get_objective(), sol2.get_objective()) - OBJECTIVE_EPSILON); + } else if (sol1.get_feasible()) { + better_than_both = offspring.get_objective() < (sol1.get_objective() - OBJECTIVE_EPSILON); + } else if (sol2.get_feasible()) { + better_than_both = offspring.get_objective() < (sol2.get_objective() - OBJECTIVE_EPSILON); + } else { + better_than_both = offspring.get_feasible(); + } + if (offspring.get_feasible() && better_than_both) { + context.settings.benchmark_info_ptr->last_improvement_after_recombination = + timer.elapsed_time(); + } +} + template std::pair, solution_t> diversity_manager_t::recombine_and_local_search(solution_t& sol1, @@ -521,9 +628,19 @@ diversity_manager_t::recombine_and_local_search(solution_t& sol1.get_feasible(), sol2.get_quality(population.weights), sol2.get_feasible()); + double best_of_parents = std::min(sol1.get_objective(), sol2.get_objective()); + bool at_least_one_parent_feasible = sol1.get_feasible() || sol2.get_feasible(); // randomly choose among 3 recombiners auto [offspring, success] = recombine(sol1, sol2); - if (!success) { return std::make_pair(solution_t(sol1), solution_t(sol2)); } + if (!success) { + // add the attempt + add_mab_reward(recombine_stats.get_last_attempt(), + std::numeric_limits::lowest(), + std::numeric_limits::lowest(), + std::numeric_limits::max(), + 0.0); + return std::make_pair(solution_t(sol1), solution_t(sol2)); + } cuopt_assert(population.test_invariant(), ""); cuopt_func_call(offspring.test_variable_bounds(false)); CUOPT_LOG_DEBUG("Recombiner offspring sol cost:feas %f : %d", @@ -531,7 +648,8 @@ diversity_manager_t::recombine_and_local_search(solution_t& offspring.get_feasible()); cuopt_assert(offspring.test_number_all_integer(), "All must be integers before LS"); bool feasibility_before = offspring.get_feasible(); - ls.run_local_search(offspring, population.weights, timer); + ls.run_local_search( + offspring, population.weights, timer, best_of_parents, at_least_one_parent_feasible); cuopt_assert(offspring.test_number_all_integer(), "All must be integers after LS"); cuopt_assert(population.test_invariant(), ""); @@ -544,14 +662,22 @@ diversity_manager_t::recombine_and_local_search(solution_t& solution_t lp_offspring(offspring); cuopt_assert(population.test_invariant(), ""); cuopt_assert(lp_offspring.test_number_all_integer(), "All must be integers before LP"); - f_t lp_run_time = offspring.get_feasible() ? 3. : 1.; + f_t lp_run_time = offspring.get_feasible() ? diversity_config_t::lp_run_time_if_feasible + : diversity_config_t::lp_run_time_if_infeasible; lp_run_time = std::min(lp_run_time, timer.remaining_time()); + relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = lp_run_time; + lp_settings.tolerance = context.settings.tolerances.absolute_tolerance; + lp_settings.return_first_feasible = false; + lp_settings.save_state = true; + lp_settings.per_constraint_residual = true; run_lp_with_vars_fixed(*lp_offspring.problem_ptr, lp_offspring, lp_offspring.problem_ptr->integer_indices, - context.settings.get_tolerances(), - context.lp_state, - lp_run_time); + lp_settings, + &ls.constraint_prop.bounds_update, + true /* check fixed assignment is feasible */, + true /* use integer fixed problem */); cuopt_assert(population.test_invariant(), ""); cuopt_assert(lp_offspring.test_number_all_integer(), "All must be integers after LP"); f_t lp_qual = lp_offspring.get_quality(population.weights); @@ -559,6 +685,16 @@ diversity_manager_t::recombine_and_local_search(solution_t& f_t offspring_qual = std::min(offspring.get_quality(population.weights), lp_qual); recombine_stats.update_improve_stats( offspring_qual, sol1.get_quality(population.weights), sol2.get_quality(population.weights)); + add_mab_reward( + recombine_stats.get_last_attempt(), + std::min(sol1.get_quality(population.weights), sol2.get_quality(population.weights)), + population.best_feasible().get_quality(population.weights), + offspring_qual, + recombine_stats.get_last_recombiner_time()); + if (context.settings.benchmark_info_ptr != nullptr) { + check_better_than_both(offspring, sol1, sol2); + check_better_than_both(lp_offspring, sol1, sol2); + } return std::make_pair(std::move(offspring), std::move(lp_offspring)); } @@ -566,26 +702,196 @@ template std::pair, bool> diversity_manager_t::recombine( solution_t& a, solution_t& b) { - i_t recombiner = std::uniform_int_distribution(0, recombiner_enum_t::SIZE - 1)(rng); + i_t recombiner; + if (run_only_ls_recombiner) { + recombiner = recombiner_enum_t::LINE_SEGMENT; + } else if (run_only_bp_recombiner) { + recombiner = recombiner_enum_t::BOUND_PROP; + } else if (run_only_fp_recombiner) { + recombiner = recombiner_enum_t::FP; + } else { + recombiner = select_mab_recombiner(); + } recombine_stats.add_attempt((recombiner_enum_t)recombiner); + recombine_stats.start_recombiner_time(); if (recombiner == recombiner_enum_t::BOUND_PROP) { - CUOPT_LOG_DEBUG("Running bound_prop recombiner"); - auto [sol, success] = bound_prop_recombiner.recombine(a, b); + auto [sol, success] = bound_prop_recombiner.recombine(a, b, population.weights); + recombine_stats.stop_recombiner_time(); if (success) { recombine_stats.add_success(); } return std::make_pair(sol, success); } else if (recombiner == recombiner_enum_t::FP) { - CUOPT_LOG_DEBUG("Running fp recombiner"); - auto [sol, success] = fp_recombiner.recombine(a, b); + auto [sol, success] = fp_recombiner.recombine(a, b, population.weights); + recombine_stats.stop_recombiner_time(); if (success) { recombine_stats.add_success(); } return std::make_pair(sol, success); } else { - CUOPT_LOG_DEBUG("Running line segment recombiner"); auto [sol, success] = line_segment_recombiner.recombine(a, b, population.weights); + recombine_stats.stop_recombiner_time(); if (success) { recombine_stats.add_success(); } return std::make_pair(sol, success); } } +template +recombiner_enum_t diversity_manager_t::select_mab_recombiner() +{ + if (recombiner_enum_t::SIZE == 0) { + CUOPT_LOG_ERROR("select_mab_recombiner called with no recombiners defined."); + cuopt_expects(false, error_type_t::RuntimeError, "No recombiners available to select in MAB."); + } + + mab_total_steps_++; + + // Phase 1: Initial exploration - ensure each arm is tried at least once + for (int i = 0; i < static_cast(recombiner_enum_t::SIZE); ++i) { + if (mab_arm_stats_[i].num_pulls == 0) { + CUOPT_LOG_DEBUG("MAB Initial Pull: Arm " + std::to_string(i)); + return static_cast(i); + } + } + + if (use_ucb_) { + // Phase 2: UCB Action Selection + return select_ucb_arm(); + } else { + // Fallback to epsilon-greedy if desired + return select_epsilon_greedy_arm(); + } +} + +// UCB arm selection with confidence bounds +template +recombiner_enum_t diversity_manager_t::select_ucb_arm() +{ + double max_ucb_value = -std::numeric_limits::infinity(); + std::vector best_arms; + + for (int i = 0; i < static_cast(recombiner_enum_t::SIZE); ++i) { + // Calculate UCB value: Q(a) + sqrt(2*ln(t)/N(a)) + double confidence_bound = std::sqrt((2.0 * std::log(mab_total_steps_)) / + static_cast(mab_arm_stats_[i].num_pulls)); + double ucb_value = mab_arm_stats_[i].q_value + confidence_bound; + + CUOPT_LOG_DEBUG( + "MAB UCB: Arm " + std::to_string(i) + ", Q=" + std::to_string(mab_arm_stats_[i].q_value) + + ", CB=" + std::to_string(confidence_bound) + ", UCB=" + std::to_string(ucb_value)); + + constexpr double tolerance = 1e-9; + if (ucb_value > max_ucb_value + tolerance) { + max_ucb_value = ucb_value; + best_arms.clear(); + best_arms.push_back(static_cast(i)); + } else if (std::abs(ucb_value - max_ucb_value) < tolerance) { + best_arms.push_back(static_cast(i)); + } + } + + if (!best_arms.empty()) { + std::uniform_int_distribution dist_tie(0, best_arms.size() - 1); + recombiner_enum_t chosen_arm = best_arms[dist_tie(mab_rng_)]; + CUOPT_LOG_DEBUG("MAB UCB Selected: Arm " + std::to_string(static_cast(chosen_arm)) + + " (UCB Value: " + std::to_string(max_ucb_value) + ")"); + return chosen_arm; + } else { + CUOPT_LOG_ERROR("MAB UCB: No best arm found, falling back to random."); + std::uniform_int_distribution dist_arm(0, recombiner_enum_t::SIZE - 1); + return static_cast(dist_arm(mab_rng_)); + } +} + +// Fallback epsilon-greedy method (preserved for compatibility) +template +recombiner_enum_t diversity_manager_t::select_epsilon_greedy_arm() +{ + std::uniform_real_distribution dist_epsilon(0.0, 1.0); + if (dist_epsilon(mab_rng_) < mab_epsilon_) { + // Explore: Choose a random arm + std::uniform_int_distribution dist_arm(0, recombiner_enum_t::SIZE - 1); + recombiner_enum_t random_arm = static_cast(dist_arm(mab_rng_)); + CUOPT_LOG_DEBUG("MAB Explore: Arm " + std::to_string(static_cast(random_arm))); + return random_arm; + } else { + // Exploit: Choose the arm with the highest Q value + double max_q_value = -std::numeric_limits::infinity(); + std::vector best_arms; + + for (int i = 0; i < static_cast(recombiner_enum_t::SIZE); ++i) { + constexpr double tolerance = 1e-9; + if (mab_arm_stats_[i].q_value > max_q_value + tolerance) { + max_q_value = mab_arm_stats_[i].q_value; + best_arms.clear(); + best_arms.push_back(static_cast(i)); + } else if (std::abs(mab_arm_stats_[i].q_value - max_q_value) < tolerance) { + best_arms.push_back(static_cast(i)); + } + } + + if (!best_arms.empty()) { + std::uniform_int_distribution dist_tie(0, best_arms.size() - 1); + recombiner_enum_t chosen_arm = best_arms[dist_tie(mab_rng_)]; + CUOPT_LOG_DEBUG("MAB Exploit: Arm " + std::to_string(static_cast(chosen_arm)) + + " (Q Value: " + std::to_string(max_q_value) + ")"); + return chosen_arm; + } + } + + // Fallback + std::uniform_int_distribution dist_arm(0, recombiner_enum_t::SIZE - 1); + return static_cast(dist_arm(mab_rng_)); +} + +template +void diversity_manager_t::add_mab_reward(recombiner_enum_t recombiner_id, + double best_of_parents_quality, + double best_feasible_quality, + double offspring_quality, + double recombination_time_in_miliseconds) +{ + int id_val = static_cast(recombiner_id); + double epsilon = std::max(1e-6, 1e-4 * std::abs(best_feasible_quality)); + bool is_better_than_best_feasible = offspring_quality + epsilon < best_feasible_quality; + bool is_better_than_best_of_parents = offspring_quality + epsilon < best_of_parents_quality; + if (id_val >= 0 && id_val < static_cast(mab_arm_stats_.size())) { + // Calculate reward based on your existing logic + double reward = 0.0; + if (is_better_than_best_feasible) { + reward = 8.0; + } else if (is_better_than_best_of_parents) { + double factor; + if (std::abs(offspring_quality - best_feasible_quality) / + (std::abs(best_feasible_quality) + 1.0) > + 1.0) { + factor = 0.; + } else if (std::abs(offspring_quality - best_feasible_quality) / + (std::abs(best_feasible_quality) + 1.0) > + 0.2) { + factor = 0.; + } else { + factor = 1.; + } + reward = factor * (std::max(0.1, 4.0 - (recombination_time_in_miliseconds / 2000))); + } + + // Update statistics + mab_arm_stats_[id_val].num_pulls++; + mab_arm_stats_[id_val].last_reward = reward; + + // Exponential recency-weighted average update: Q_new = Q_old + α(R - Q_old) + double prediction_error = reward - mab_arm_stats_[id_val].q_value; + mab_arm_stats_[id_val].q_value += mab_alpha_ * prediction_error; + + CUOPT_LOG_DEBUG( + "MAB Reward Update: Arm " + std::to_string(id_val) + ", Reward: " + std::to_string(reward) + + ", is_better_than_best_of_parents: " + (is_better_than_best_of_parents ? "Yes" : "No") + + ", Better than best: " + (is_better_than_best_feasible ? "Yes" : "No") + + ", Pulls: " + std::to_string(mab_arm_stats_[id_val].num_pulls) + + ", Q Value: " + std::to_string(mab_arm_stats_[id_val].q_value)); + } else { + CUOPT_LOG_ERROR("MAB: Attempted to add reward for invalid recombiner_id: " + + std::to_string(id_val)); + } +} + #if MIP_INSTANTIATE_FLOAT template class diversity_manager_t; #endif diff --git a/cpp/src/mip/diversity/diversity_manager.cuh b/cpp/src/mip/diversity/diversity_manager.cuh index 4a2da6bbe..c45b83107 100644 --- a/cpp/src/mip/diversity/diversity_manager.cuh +++ b/cpp/src/mip/diversity/diversity_manager.cuh @@ -43,6 +43,7 @@ class diversity_manager_t { solution_t generate_solution(f_t time_limit, bool random_start = true); // generates initial solutions void generate_initial_solutions(); + void run_fj_alone(solution_t& solution); // main loop of diversity improvements void main_loop(); // randomly chooses a recombiner and returns the offspring @@ -54,7 +55,7 @@ class diversity_manager_t { void average_fj_weights(i_t i); void diversity_step(); std::vector> generate_more_solutions(); - void add_user_given_solution(std::vector>& initial_sol_vector); + void add_user_given_solutions(std::vector>& initial_sol_vector); population_t* get_population_pointer() { return &population; } void recombine_and_ls_with_all(std::vector>& solutions); void recombine_and_ls_with_all(solution_t& solution); @@ -63,6 +64,9 @@ class diversity_manager_t { void set_new_user_bound(f_t new_user_bound); void generate_quick_feasible_solution(); bool check_b_b_preemption(); + void check_better_than_both(solution_t& offspring, + solution_t& sol1, + solution_t& sol2); mip_solver_context_t& context; problem_t* problem_ptr; @@ -78,6 +82,32 @@ class diversity_manager_t { i_t current_step{0}; solver_stats_t& stats; std::vector> initial_sol_vector; + + // Enhanced statistics structure for UCB with exponential recency weighting + struct mab_arm_stats_t { + int num_pulls = 0; // Number of times this arm was selected + double q_value = 0.5; // Exponential recency-weighted average estimate + double last_reward = 0.0; // Last reward received (for debugging) + }; + std::vector mab_arm_stats_; + double mab_epsilon_ = 0.15; // Probability of exploration in Epsilon-Greedy. + std::mt19937 mab_rng_; // RNG dedicated to MAB decisions. + double mab_alpha_ = 0.05; // Step size for exponential recency weighting + int mab_total_steps_ = 0; // Total number of action selections (for UCB) + bool use_ucb_ = true; // Flag to enable UCB vs epsilon-greedy + + // --- MAB Helper Methods --- + recombiner_enum_t select_mab_recombiner(); + void add_mab_reward(recombiner_enum_t recombiner_id, + double best_of_parents_quality, + double best_feasible_quality, + double offspring_quality, + double recombination_time_in_miliseconds); + recombiner_enum_t select_ucb_arm(); + recombiner_enum_t select_epsilon_greedy_arm(); + bool run_only_ls_recombiner{false}; + bool run_only_bp_recombiner{false}; + bool run_only_fp_recombiner{false}; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/diversity/population.cu b/cpp/src/mip/diversity/population.cu index 0b4001772..5ed21f9c4 100644 --- a/cpp/src/mip/diversity/population.cu +++ b/cpp/src/mip/diversity/population.cu @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -27,10 +28,12 @@ namespace cuopt::linear_programming::detail { -constexpr double weight_increase_ratio = 2.; -constexpr double weight_decrease_ratio = 0.9; -constexpr double max_infeasibility_weight = 10000000.; -constexpr double min_infeasibility_weight = 1.; +constexpr double weight_increase_ratio = 2.; +constexpr double weight_decrease_ratio = 0.9; +constexpr double max_infeasibility_weight = 1e12; +constexpr double min_infeasibility_weight = 1.; +constexpr double infeasibility_balance_ratio = 1.1; +constexpr double halving_skip_ratio = 0.75; template population_t::population_t(std::string const& name_, @@ -46,10 +49,10 @@ population_t::population_t(std::string const& name_, infeasibility_importance(infeasibility_weight_), weights(0, context.problem_ptr->handle_ptr), rng(cuopt::seed_generator::get_seed()), - early_exit_primal_generation(false) + early_exit_primal_generation(false), + timer(0) { - best_feasible_objective = - problem_ptr->maximize ? -std::numeric_limits::max() : std::numeric_limits::max(); + best_feasible_objective = std::numeric_limits::max(); } template @@ -66,7 +69,6 @@ void population_t::initialize_population() { var_threshold = std::max(problem_ptr->n_variables - var_threshold, (problem_ptr->n_variables / 10) * 8); - initial_threshold_ratio = (f_t)var_threshold / problem_ptr->n_variables; solutions.reserve(max_solutions); indices.reserve(max_solutions); // indices[0] always points to solutions[0] - a special place for feasible solution @@ -85,7 +87,7 @@ std::pair, solution_t> population_t::ge { raft::common::nvtx::range fun_scope("get_two_random"); cuopt_assert(indices.size() > 2, "There should be enough solutions"); - size_t add = (size_t)(!solutions[0].first || solutions[indices[1].first].second.get_feasible()); + size_t add = (size_t)(!solutions[0].first); size_t i = add + std::uniform_int_distribution(0, (indices.size() - 2))(rng); size_t j = add + std::uniform_int_distribution(0, (indices.size() - 3))(rng); if (tournament) { @@ -97,6 +99,19 @@ std::pair, solution_t> population_t::ge if (j >= i) j++; auto first_solution = solutions[indices[i].first].second; auto second_solution = solutions[indices[j].first].second; + // if best feasible and best are the same, take the second index instead of best + if (i == 0 && j == 1) { + bool same = + check_integer_equal_on_indices(first_solution.problem_ptr->integer_indices, + first_solution.assignment, + second_solution.assignment, + first_solution.problem_ptr->tolerances.integrality_tolerance, + first_solution.handle_ptr); + if (same) { + auto new_sol = solutions[indices[2].first].second; + second_solution = std::move(new_sol); + } + } cuopt_assert(test_invariant(), "Population invariant doesn't hold"); return std::make_pair(std::move(first_solution), std::move(second_solution)); } @@ -157,16 +172,23 @@ std::vector> population_t::get_external_solutions return return_vector; } +template +bool population_t::is_better_than_best_feasible(solution_t& sol) +{ + bool obj_better = sol.get_objective() < best_feasible_objective; + return obj_better && sol.get_feasible(); +} + template void population_t::run_solution_callbacks(solution_t& sol) { - bool better_solution_found = problem_ptr->maximize - ? sol.get_user_objective() > best_feasible_objective - : sol.get_user_objective() < best_feasible_objective; + bool better_solution_found = is_better_than_best_feasible(sol); auto user_callbacks = context.settings.get_mip_callbacks(); - if (better_solution_found && sol.get_feasible()) { + if (better_solution_found) { + if (context.settings.benchmark_info_ptr != nullptr) { + context.settings.benchmark_info_ptr->last_improvement_of_best_feasible = timer.elapsed_time(); + } CUOPT_LOG_DEBUG("Population: Found new best solution %g", sol.get_user_objective()); - best_feasible_objective = sol.get_user_objective(); if (problem_ptr->branch_and_bound_callback != nullptr) { problem_ptr->branch_and_bound_callback(sol.get_host_assignment()); } @@ -199,6 +221,11 @@ void population_t::run_solution_callbacks(solution_t& sol) get_sol_callback->get_solution(temp_sol.assignment.data(), user_objective_vec.data()); } } + // save the best objective here, because we might not have been able to return the solution to + // the user because of the unscaling that causes infeasibility. + // This prevents an issue of repaired, or a fully feasible solution being reported in the call + // back in next run. + best_feasible_objective = sol.get_objective(); } for (auto callback : user_callbacks) { @@ -245,6 +272,38 @@ void population_t::run_solution_callbacks(solution_t& sol) } } +template +void population_t::adjust_weights_according_to_best_feasible() +{ + // check if the best in population still the best feasible + if (!best().get_feasible()) { + CUOPT_LOG_DEBUG("Best solution is infeasible, adjusting weights"); + // if not the case, adjust the weights such that the best is feasible + f_t weighted_violation_of_best = best().get_quality(weights) - best().get_objective(); + CUOPT_LOG_DEBUG("weighted_violation_of_best %f quality %f objective %f\n", + weighted_violation_of_best, + best().get_quality(weights), + best().get_objective()); + cuopt_assert(weighted_violation_of_best > 1e-10, "Weighted violation of best is not positive"); + // fixme + weighted_violation_of_best = max(weighted_violation_of_best, 1e-10); + f_t quality_difference = best_feasible().get_quality(weights) - best().get_quality(weights); + CUOPT_LOG_DEBUG("quality_difference %f best_feasible_quality %f best_quality %f", + quality_difference, + best_feasible().get_quality(weights), + best().get_quality(weights)); + if (quality_difference < 1e-10) { return; } + // make the current best infeasible 10% worse than feasible + f_t increase_ratio = + (quality_difference * infeasibility_balance_ratio) / weighted_violation_of_best; + infeasibility_importance *= (1 + increase_ratio); + infeasibility_importance = min(max_infeasibility_weight, infeasibility_importance); + normalize_weights(); + update_qualities(); + cuopt_assert(test_invariant(), "Population invariant doesn't hold"); + } +} + template i_t population_t::add_solution(solution_t&& sol) { @@ -408,17 +467,6 @@ void population_t::compute_new_weights() best_sol.handle_ptr->sync_stream(); } -template -void population_t::adjust_threshold(cuopt::timer_t timer) -{ - const double max_diversity_threshold = 0.99; - double time_ratio = - (timer.elapsed_time() - diversity_start_time) / (timer.get_time_limit() - diversity_start_time); - f_t threshold_ratio = - initial_threshold_ratio + time_ratio * (max_diversity_threshold - initial_threshold_ratio); - var_threshold = threshold_ratio * problem_ptr->n_variables; -} - template void population_t::update_qualities() { @@ -439,9 +487,6 @@ void population_t::update_weights() { raft::common::nvtx::range fun_scope("adjust_weight_changes"); CUOPT_LOG_DEBUG("Changing the weights"); - // TODO activate this if we have a reserve and a diverse initial population - // by adding new solutions at every diversity step, it doesn't make sense to add - // adjust_threshold(timer); compute_new_weights(); normalize_weights(); update_qualities(); @@ -550,7 +595,8 @@ template void population_t::halve_the_population() { raft::common::nvtx::range fun_scope("halve_the_population"); - if (current_size() <= max_solutions / 2) { return; } + // try 3/4 here + if (current_size() <= (max_solutions * halving_skip_ratio)) { return; } CUOPT_LOG_DEBUG("Halving the population, current size: %lu", current_size()); // put population into a vector auto sol_vec = population_to_vector(); @@ -571,7 +617,7 @@ void population_t::halve_the_population() clear_except_best_feasible(); var_threshold = std::min( max_var_threshold, - std::min((size_t)(var_threshold * 0.97), (size_t)(0.995 * problem_ptr->n_integer_vars))); + std::min((size_t)(var_threshold * 1.02), (size_t)(0.995 * problem_ptr->n_integer_vars))); for (auto& sol : sol_vec) { add_solution(solution_t(sol)); } @@ -591,6 +637,23 @@ size_t population_t::find_free_solution_index() return std::numeric_limits::max(); } +template +void population_t::start_threshold_adjustment() +{ + population_start_time = timer.elapsed_time(); + initial_threshold = var_threshold; +} + +template +void population_t::adjust_threshold(cuopt::timer_t timer) +{ + double time_ratio = (timer.elapsed_time() - population_start_time) / + (timer.get_time_limit() - population_start_time); + var_threshold = + initial_threshold + + time_ratio * (get_max_var_threshold(problem_ptr->n_integer_vars) - initial_threshold); +} + template void population_t::find_diversity(std::vector>& initial_sol_vector, bool avg) @@ -682,11 +745,12 @@ void population_t::print() if (index.first == 0 && solutions[0].first) { CUOPT_LOG_DEBUG(" Best feasible: %f", solutions[index.first].second.get_user_objective()); } - CUOPT_LOG_DEBUG("%d : %f\t%f\t%f", + CUOPT_LOG_DEBUG("%d : %f\t%f\t%f\t%d", i, index.second, solutions[index.first].second.get_total_excess(), - solutions[index.first].second.get_user_objective()); + solutions[index.first].second.get_user_objective(), + solutions[index.first].second.get_feasible()); i++; } CUOPT_LOG_DEBUG(" -------------- "); diff --git a/cpp/src/mip/diversity/population.cuh b/cpp/src/mip/diversity/population.cuh index 0d103b93e..0f0176341 100644 --- a/cpp/src/mip/diversity/population.cuh +++ b/cpp/src/mip/diversity/population.cuh @@ -63,6 +63,7 @@ class population_t { solution_t& solution_at_index(i_t idx) { return solutions[indices[idx].first].second; } // initializes the population lazily. after presolve and var removals void initialize_population(); + bool is_better_than_best_feasible(solution_t& sol); void allocate_solutions(); @@ -141,6 +142,10 @@ class population_t { void run_solution_callbacks(solution_t& sol); + void adjust_weights_according_to_best_feasible(); + + void start_threshold_adjustment(); + // does some consistency tests bool test_invariant(); @@ -150,8 +155,8 @@ class population_t { mip_solver_context_t& context; problem_t* problem_ptr; i_t var_threshold; - f_t initial_threshold_ratio; - f_t diversity_start_time; + i_t initial_threshold; + double population_start_time; // the normalization target for the infeasibility // this is used to cover the importance of the weights f_t infeasibility_importance = 100.; @@ -166,6 +171,7 @@ class population_t { bool early_exit_primal_generation = false; f_t best_feasible_objective = std::numeric_limits::max(); bool preempt_heuristic_solver_ = false; + cuopt::timer_t timer; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh b/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh index be73d5dc7..f38cc5759 100644 --- a/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh +++ b/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh @@ -35,96 +35,219 @@ class bound_prop_recombiner_t : public recombiner_t { const raft::handle_t* handle_ptr) : recombiner_t(context, n_vars, handle_ptr), constraint_prop(constraint_prop_), - rng(cuopt::seed_generator::get_seed()) + rng(cuopt::seed_generator::get_seed()), + vars_to_fix(n_vars, handle_ptr->get_stream()) { } - std::pair, bool> recombine(solution_t& a, solution_t& b) + void get_probing_values_for_infeasible( + solution_t& guiding, + solution_t& other, + solution_t& offspring, + rmm::device_uvector>& probing_values, + i_t n_vars_from_other) + { + auto guiding_view = guiding.view(); + auto other_view = other.view(); + auto offspring_view = offspring.view(); + const f_t int_tol = guiding.problem_ptr->tolerances.integrality_tolerance; + // this is to give two possibilities to round in case of conflict + thrust::for_each( + guiding.handle_ptr->get_thrust_policy(), + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(guiding.problem_ptr->n_variables), + [guiding_view, other_view, probing_values = probing_values.data()] __device__(i_t idx) { + f_t guiding_val = guiding_view.assignment[idx]; + f_t other_val = other_view.assignment[idx]; + cuopt_assert(guiding_view.problem.check_variable_within_bounds(idx, guiding_val), ""); + cuopt_assert(other_view.problem.check_variable_within_bounds(idx, other_val), ""); + probing_values[idx] = thrust::make_pair(guiding_val, other_val); + }); + // populate remaining N integers randomly/average from each solution + thrust::for_each( + guiding.handle_ptr->get_thrust_policy(), + this->remaining_indices.data(), + this->remaining_indices.data() + n_vars_from_other, + [guiding_view, + other_view, + offspring_view, + int_tol, + probing_values = probing_values.data(), + seed = cuopt::seed_generator::get_seed()] __device__(i_t idx) { + f_t guiding_val = guiding_view.assignment[idx]; + f_t other_val = other_view.assignment[idx]; + cuopt_assert(guiding_view.problem.check_variable_within_bounds(idx, guiding_val), ""); + cuopt_assert(other_view.problem.check_variable_within_bounds(idx, other_val), ""); + f_t avg_val = (other_val + guiding_val) / 2; + if (guiding_view.problem.is_integer_var(idx)) { + raft::random::PCGenerator rng(seed, idx, 0); + if (rng.next_u32() % 2) { cuda::std::swap(other_val, guiding_val); } + cuopt_assert(is_integer(other_val, int_tol), "The value must be integer"); + f_t second_val = round(avg_val) == other_val ? guiding_val : round(avg_val); + probing_values[idx] = thrust::make_pair(other_val, second_val); + // assign some floating value, so that they can be rounded by bounds prop + f_t lb = guiding_view.problem.variable_lower_bounds[idx]; + f_t ub = guiding_view.problem.variable_upper_bounds[idx]; + if (integer_equal(lb, ub, int_tol)) { + cuopt_assert(false, "The var values must be different in A and B!"); + } else if (isfinite(lb)) { + offspring_view.assignment[idx] = lb + 0.1; + } else { + offspring_view.assignment[idx] = ub - 0.1; + } + } else { + // if the var is continuous, take the average + offspring_view.assignment[idx] = avg_val; + } + }); + } + + void get_probing_values_for_feasible(solution_t& guiding, + solution_t& other, + solution_t& offspring, + rmm::device_uvector>& probing_values, + i_t n_vars_from_other, + rmm::device_uvector& variable_map) + { + cuopt_assert(n_vars_from_other == offspring.problem_ptr->n_integer_vars, + "The number of vars from other should match!"); + auto guiding_view = guiding.view(); + auto other_view = other.view(); + auto offspring_view = offspring.view(); + const f_t int_tol = guiding.problem_ptr->tolerances.integrality_tolerance; + thrust::for_each( + guiding.handle_ptr->get_thrust_policy(), + thrust::make_counting_iterator(0lu), + thrust::make_counting_iterator(variable_map.size()), + [guiding_view, + other_view, + offspring_view, + int_tol, + probing_values = make_span(probing_values), + variable_map = make_span(variable_map)] __device__(size_t idx) { + f_t other_val = other_view.assignment[variable_map[idx]]; + f_t guiding_val = guiding_view.assignment[variable_map[idx]]; + cuopt_assert(other_view.problem.check_variable_within_bounds(variable_map[idx], other_val), + ""); + cuopt_assert( + guiding_view.problem.check_variable_within_bounds(variable_map[idx], guiding_val), ""); + f_t avg_val = (other_val + guiding_val) / 2; + probing_values[idx] = thrust::make_pair(guiding_val, other_val); + offspring_view.assignment[idx] = avg_val; + }); + } + + std::pair, bool> recombine(solution_t& a, + solution_t& b, + const weight_t& weights) { raft::common::nvtx::range fun_scope("bound_prop_recombiner"); - // copy the solution from A - solution_t offspring(a); + auto& guiding_solution = a.get_feasible() ? a : b; + auto& other_solution = a.get_feasible() ? b : a; + // copy the solution from guiding + solution_t offspring(guiding_solution); // find same values and populate it to offspring - this->assign_same_integer_values(a, b, offspring); - i_t remaining_variables = this->n_remaining.value(a.handle_ptr->get_stream()); - // // from the remaining integers, populate randomly. - CUOPT_LOG_DEBUG("n_vars from A/B %d remaining_variables %d", - a.problem_ptr->n_variables - remaining_variables, - remaining_variables); + i_t n_different_vars = this->assign_same_integer_values(a, b, offspring); + CUOPT_LOG_DEBUG("BP rec: Number of different variables %d MAX_VARS %d", + n_different_vars, + bp_recombiner_config_t::max_n_of_vars_from_other); + i_t n_vars_from_other = n_different_vars; + i_t fixed_from_guiding = 0; + i_t fixed_from_other = 0; + if (n_different_vars > (i_t)bp_recombiner_config_t::max_n_of_vars_from_other) { + fixed_from_guiding = n_vars_from_other - bp_recombiner_config_t::max_n_of_vars_from_other; + n_vars_from_other = bp_recombiner_config_t::max_n_of_vars_from_other; + thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()}; + thrust::shuffle(a.handle_ptr->get_thrust_policy(), + this->remaining_indices.data(), + this->remaining_indices.data() + n_different_vars, + g); + } + i_t n_vars_from_guiding = a.problem_ptr->n_integer_vars - n_vars_from_other; + CUOPT_LOG_DEBUG( + "n_vars_from_guiding %d n_vars_from_other %d", n_vars_from_guiding, n_vars_from_other); // if either all integers are from A(meaning all are common) or all integers are from B(meaning // all are different), return - if (a.problem_ptr->n_integer_vars - remaining_variables == 0 || remaining_variables == 0) { + if (n_vars_from_guiding == 0 || n_vars_from_other == 0) { + CUOPT_LOG_DEBUG("Returning false because all vars are common or different"); return std::make_pair(offspring, false); } cuopt_assert(a.problem_ptr == b.problem_ptr, "The two solutions should not refer to different problems"); - auto a_view = a.view(); - auto b_view = b.view(); - auto offspring_view = offspring.view(); - const f_t int_tol = a.problem_ptr->tolerances.integrality_tolerance; + const f_t lp_run_time_after_feasible = bp_recombiner_config_t::lp_after_bounds_prop_time_limit; + constraint_prop.max_n_failed_repair_iterations = bp_recombiner_config_t::n_repair_iterations; rmm::device_uvector> probing_values(a.problem_ptr->n_variables, a.handle_ptr->get_stream()); - // this is to give two possibilities to round in case of conflict - thrust::for_each(a.handle_ptr->get_thrust_policy(), - thrust::make_counting_iterator(0), - thrust::make_counting_iterator(a.problem_ptr->n_variables), - [a_view, b_view, probing_values = probing_values.data()] __device__(i_t idx) { - f_t a_val = a_view.assignment[idx]; - f_t b_val = b_view.assignment[idx]; - cuopt_assert(a_view.problem.check_variable_within_bounds(idx, a_val), ""); - cuopt_assert(b_view.problem.check_variable_within_bounds(idx, b_val), ""); - probing_values[idx] = thrust::make_pair(a_val, b_val); - }); - // populate remaining N integers randomly/average from each solution - thrust::for_each(a.handle_ptr->get_thrust_policy(), - this->remaining_indices.data(), - this->remaining_indices.data() + remaining_variables, - [a_view, - b_view, - offspring_view, - int_tol, - probing_values = probing_values.data(), - seed = cuopt::seed_generator::get_seed()] __device__(i_t idx) { - raft::random::PCGenerator rng(seed, idx, 0); - f_t a_val = a_view.assignment[idx]; - f_t b_val = b_view.assignment[idx]; - cuopt_assert(a_view.problem.check_variable_within_bounds(idx, a_val), ""); - cuopt_assert(b_view.problem.check_variable_within_bounds(idx, b_val), ""); - f_t avg_val = (b_val + a_val) / 2; - if (a_view.problem.is_integer_var(idx)) { - const bool rnd = rng.next_u32() % 2; - // if the var is integer, populate probing vals - f_t first_val = rnd ? b_val : a_val; - cuopt_assert(is_integer(first_val, int_tol), - "The value must be integer"); - // TODO check the rounding direction and var bounds - probing_values[idx] = thrust::make_pair(first_val, round(avg_val)); - // assign some floating value, so that they can be rounded by bounds prop - f_t lb = a_view.problem.variable_lower_bounds[idx]; - f_t ub = a_view.problem.variable_upper_bounds[idx]; - if (integer_equal(lb, ub, int_tol)) { - cuopt_assert(false, "The var values must be different in A and B!"); - } else if (isfinite(lb)) { - offspring_view.assignment[idx] = lb + 0.1; - } else { - offspring_view.assignment[idx] = ub - 0.1; - } - } else { - // if the var is continuous, take the average - offspring_view.assignment[idx] = avg_val; - } - }); - const f_t lp_run_time_after_feasible = 2.; - timer_t timer(2.); - auto h_probing_values = host_copy(probing_values); - constraint_prop.apply_round(offspring, lp_run_time_after_feasible, timer, h_probing_values); + probing_config_t probing_config(a.problem_ptr->n_variables, a.handle_ptr); + if (guiding_solution.get_feasible()) { + this->compute_vars_to_fix(offspring, vars_to_fix, n_vars_from_other, n_vars_from_guiding); + auto [fixed_problem, fixed_assignment, variable_map] = offspring.fix_variables(vars_to_fix); + timer_t timer(bp_recombiner_config_t::bounds_prop_time_limit); + rmm::device_uvector old_assignment(offspring.assignment, + offspring.handle_ptr->get_stream()); + offspring.handle_ptr->sync_stream(); + offspring.assignment = std::move(fixed_assignment); + offspring.problem_ptr = &fixed_problem; + cuopt_func_call(offspring.test_variable_bounds(false)); + get_probing_values_for_feasible(guiding_solution, + other_solution, + offspring, + probing_values, + n_vars_from_other, + variable_map); + probing_config.probing_values = host_copy(probing_values); + probing_config.n_of_fixed_from_first = fixed_from_guiding; + probing_config.n_of_fixed_from_second = fixed_from_other; + probing_config.use_balanced_probing = true; + constraint_prop.single_rounding_only = true; + constraint_prop.apply_round(offspring, lp_run_time_after_feasible, timer, probing_config); + constraint_prop.single_rounding_only = false; + cuopt_func_call(bool feasible_after_bounds_prop = offspring.get_feasible()); + offspring.handle_ptr->sync_stream(); + offspring.problem_ptr = a.problem_ptr; + fixed_assignment = std::move(offspring.assignment); + offspring.assignment = std::move(old_assignment); + 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!"); + a.handle_ptr->sync_stream(); + } else { + timer_t timer(bp_recombiner_config_t::bounds_prop_time_limit); + get_probing_values_for_infeasible( + guiding_solution, other_solution, offspring, probing_values, n_vars_from_other); + probing_config.probing_values = host_copy(probing_values); + constraint_prop.apply_round(offspring, lp_run_time_after_feasible, timer, probing_config); + } + constraint_prop.max_n_failed_repair_iterations = 1; cuopt_func_call(offspring.test_number_all_integer()); - offspring.compute_feasibility(); - bool same_as_parents = this->check_if_offspring_is_same_as_parents(offspring, a, b); + bool better_cost_than_parents = + offspring.get_quality(weights) < + std::min(other_solution.get_quality(weights), guiding_solution.get_quality(weights)); + bool better_feasibility_than_parents = offspring.get_feasible() && + !other_solution.get_feasible() && + !guiding_solution.get_feasible(); + + bool same_as_parents = + this->check_if_offspring_is_same_as_parents(offspring, guiding_solution, other_solution); + // adjust the max_n_of_vars_from_other + if (n_different_vars > (i_t)bp_recombiner_config_t::max_n_of_vars_from_other) { + if (same_as_parents) { + bp_recombiner_config_t::increase_max_n_of_vars_from_other(); + } else { + bp_recombiner_config_t::decrease_max_n_of_vars_from_other(); + } + } + if (better_cost_than_parents || better_feasibility_than_parents) { + CUOPT_LOG_DEBUG("Offspring is feasible or better than both parents"); + return std::make_pair(offspring, true); + } return std::make_pair(offspring, !same_as_parents); } + rmm::device_uvector vars_to_fix; constraint_prop_t& constraint_prop; thrust::default_random_engine rng; }; diff --git a/cpp/src/mip/diversity/recombiners/fp_recombiner.cuh b/cpp/src/mip/diversity/recombiners/fp_recombiner.cuh index aa6206ff0..5597e8e84 100644 --- a/cpp/src/mip/diversity/recombiners/fp_recombiner.cuh +++ b/cpp/src/mip/diversity/recombiners/fp_recombiner.cuh @@ -26,8 +26,6 @@ #include #include -#include -#include namespace cuopt::linear_programming::detail { @@ -44,70 +42,57 @@ class fp_recombiner_t : public recombiner_t { { } - std::pair, bool> recombine(solution_t& a, solution_t& b) + std::pair, bool> recombine(solution_t& a, + solution_t& b, + const weight_t& weights) { raft::common::nvtx::range fun_scope("FP recombiner"); - + auto& guiding_solution = a.get_feasible() ? a : b; + auto& other_solution = a.get_feasible() ? b : a; // copy the solution from A - solution_t offspring(a); + solution_t offspring(guiding_solution); // find same values and populate it to offspring - this->assign_same_integer_values(a, b, offspring); - const i_t n_remaining_vars = this->n_remaining.value(a.handle_ptr->get_stream()); - // partition to get only integer uncommon - auto iter = thrust::stable_partition(a.handle_ptr->get_thrust_policy(), - this->remaining_indices.begin(), - this->remaining_indices.begin() + n_remaining_vars, - [a_view = a.view()] __device__(i_t idx) { - if (a_view.problem.is_integer_var(idx)) { return true; } - return false; - }); - i_t h_n_remaining_integers = iter - this->remaining_indices.begin(); - i_t n_common_integers = a.problem_ptr->n_integer_vars - h_n_remaining_integers; - if (h_n_remaining_integers == 0 || n_common_integers == 0) { - CUOPT_LOG_DEBUG("All integers are common or different in FP recombiner, returning A!"); + i_t n_different_vars = + this->assign_same_integer_values(guiding_solution, other_solution, offspring); + CUOPT_LOG_DEBUG("FP rec: Number of different variables %d MAX_VARS %d", + n_different_vars, + fp_recombiner_config_t::max_n_of_vars_from_other); + i_t n_vars_from_other = n_different_vars; + if (n_vars_from_other > (i_t)fp_recombiner_config_t::max_n_of_vars_from_other) { + n_vars_from_other = fp_recombiner_config_t::max_n_of_vars_from_other; + thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()}; + thrust::shuffle(a.handle_ptr->get_thrust_policy(), + this->remaining_indices.data(), + this->remaining_indices.data() + n_different_vars, + g); + } + i_t n_vars_from_guiding = a.problem_ptr->n_integer_vars - n_vars_from_other; + if (n_vars_from_other == 0 || n_vars_from_guiding == 0) { + CUOPT_LOG_DEBUG("Returning false because all vars are common or different"); return std::make_pair(offspring, false); } - CUOPT_LOG_DEBUG("n_integer_vars from A/B %d n_integer_vars from B %d", - n_common_integers, - h_n_remaining_integers); - vars_to_fix.resize(n_common_integers, a.handle_ptr->get_stream()); - // set difference needs two sorted arrays - thrust::sort(a.handle_ptr->get_thrust_policy(), - this->remaining_indices.data(), - this->remaining_indices.data() + h_n_remaining_integers); - cuopt_assert((thrust::is_sorted(a.handle_ptr->get_thrust_policy(), - a.problem_ptr->integer_indices.begin(), - a.problem_ptr->integer_indices.end())), - "vars_to_fix should be sorted!"); - // get the variables to fix (common variables) - iter = thrust::set_difference(a.handle_ptr->get_thrust_policy(), - a.problem_ptr->integer_indices.begin(), - a.problem_ptr->integer_indices.end(), - this->remaining_indices.data(), - this->remaining_indices.data() + h_n_remaining_integers, - vars_to_fix.begin()); - cuopt_assert(iter - vars_to_fix.begin() == n_common_integers, "The size should match!"); - cuopt_assert((thrust::is_sorted(a.handle_ptr->get_thrust_policy(), - vars_to_fix.data(), - vars_to_fix.data() + n_common_integers)), - "vars_to_fix should be sorted!"); - const f_t tolerance = 1e-2; - const f_t lp_time = 0.5; + CUOPT_LOG_DEBUG( + "n_vars_from_guiding %d n_vars_from_other %d", n_vars_from_guiding, n_vars_from_other); + this->compute_vars_to_fix(offspring, vars_to_fix, n_vars_from_other, n_vars_from_guiding); auto [fixed_problem, fixed_assignment, variable_map] = offspring.fix_variables(vars_to_fix); fixed_problem.check_problem_representation(true); - const bool check_feas = false; - const bool return_first_feasible = false; - const bool save_state = false; - // every sub problem is different,so it is very hard to find a valid initial solution - lp_state_t lp_state = this->context.lp_state; - auto solver_response = get_relaxed_lp_solution(fixed_problem, - fixed_assignment, - lp_state, - tolerance, - lp_time, - check_feas, - return_first_feasible, - save_state); + if (!guiding_solution.get_feasible() && !other_solution.get_feasible()) { + relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = fp_recombiner_config_t::infeasibility_detection_time_limit; + lp_settings.tolerance = fixed_problem.tolerances.absolute_tolerance; + lp_settings.return_first_feasible = true; + lp_settings.save_state = true; + lp_settings.check_infeasibility = true; + // run lp with infeasibility detection on + auto lp_response = + get_relaxed_lp_solution(fixed_problem, fixed_assignment, offspring.lp_state, lp_settings); + if (lp_response.get_termination_status() == pdlp_termination_status_t::PrimalInfeasible || + lp_response.get_termination_status() == pdlp_termination_status_t::DualInfeasible || + lp_response.get_termination_status() == pdlp_termination_status_t::TimeLimit) { + CUOPT_LOG_DEBUG("FP recombiner failed because LP found infeasible!"); + return std::make_pair(offspring, false); + } + } // brute force rounding threshold is 8 const bool run_fp = fixed_problem.n_integer_vars > 8; if (run_fp) { @@ -117,15 +102,16 @@ class fp_recombiner_t : public recombiner_t { offspring.handle_ptr->get_stream()); offspring.handle_ptr->sync_stream(); offspring.assignment = std::move(fixed_assignment); - offspring.lp_state = std::move(lp_state); cuopt_func_call(offspring.test_variable_bounds(false)); - timer_t timer((f_t)2.); + timer_t timer(fp_recombiner_config_t::fp_time_limit); fp.timer = timer; fp.cycle_queue.reset(offspring); fp.reset(); fp.resize_vectors(*offspring.problem_ptr, offspring.handle_ptr); - bool is_feasible = fp.run_single_fp_descent(offspring); - if (is_feasible) { CUOPT_LOG_DEBUG("FP after recombiner, found feasible!"); } + fp.config.alpha = fp_recombiner_config_t::alpha; + fp.config.alpha_decrease_factor = fp_recombiner_config_t::alpha_decrease_factor; + bool is_feasible = fp.run_single_fp_descent(offspring); + if (is_feasible) { CUOPT_LOG_DEBUG("FP recombiner found feasible!"); } CUOPT_LOG_DEBUG("FP completed after recombiner!"); offspring.handle_ptr->sync_stream(); offspring.problem_ptr = orig_problem_ptr; @@ -139,6 +125,24 @@ class fp_recombiner_t : public recombiner_t { cuopt_assert(offspring.test_number_all_integer(), "All must be integers after offspring"); offspring.compute_feasibility(); bool same_as_parents = this->check_if_offspring_is_same_as_parents(offspring, a, b); + // adjust the max_n_of_vars_from_other + if (n_different_vars > (i_t)fp_recombiner_config_t::max_n_of_vars_from_other) { + if (same_as_parents) { + fp_recombiner_config_t::increase_max_n_of_vars_from_other(); + } else { + fp_recombiner_config_t::decrease_max_n_of_vars_from_other(); + } + } + bool better_cost_than_parents = + offspring.get_quality(weights) < + std::min(other_solution.get_quality(weights), guiding_solution.get_quality(weights)); + bool better_feasibility_than_parents = offspring.get_feasible() && + !other_solution.get_feasible() && + !guiding_solution.get_feasible(); + if (better_cost_than_parents || better_feasibility_than_parents) { + CUOPT_LOG_DEBUG("Offspring is feasible or better than both parents"); + return std::make_pair(offspring, true); + } return std::make_pair(offspring, !same_as_parents); } rmm::device_uvector vars_to_fix; diff --git a/cpp/src/mip/diversity/recombiners/line_segment_recombiner.cuh b/cpp/src/mip/diversity/recombiners/line_segment_recombiner.cuh index 78c9afc34..cfde6d03c 100644 --- a/cpp/src/mip/diversity/recombiners/line_segment_recombiner.cuh +++ b/cpp/src/mip/diversity/recombiners/line_segment_recombiner.cuh @@ -36,28 +36,99 @@ class line_segment_recombiner_t : public recombiner_t { : recombiner_t(context, n_vars, handle_ptr), line_segment_search(line_segment_search_) { } + + rmm::device_uvector generate_delta_vector(solution_t& guiding_solution, + solution_t& other_solution, + solution_t& offspring, + i_t n_points_to_search, + i_t remaining_variables) + { + CUOPT_LOG_DEBUG("LS rec: Number of different variables %d MAX_VARS %d", + remaining_variables, + ls_recombiner_config_t::max_n_of_vars_from_other); + i_t n_vars_from_other = remaining_variables; + if (n_vars_from_other > (i_t)ls_recombiner_config_t::max_n_of_vars_from_other) { + n_vars_from_other = ls_recombiner_config_t::max_n_of_vars_from_other; + thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()}; + thrust::shuffle(guiding_solution.handle_ptr->get_thrust_policy(), + this->remaining_indices.data(), + this->remaining_indices.data() + remaining_variables, + g); + } + i_t n_vars_from_guiding = guiding_solution.problem_ptr->n_integer_vars - n_vars_from_other; + rmm::device_uvector delta_vector(offspring.problem_ptr->n_variables, + offspring.handle_ptr->get_stream()); + thrust::fill( + offspring.handle_ptr->get_thrust_policy(), delta_vector.begin(), delta_vector.end(), (f_t)0); + // generate delta vector only for the ones we want to search + thrust::for_each(offspring.handle_ptr->get_thrust_policy(), + this->remaining_indices.data(), + this->remaining_indices.data() + n_vars_from_other, + [guiding_solution = guiding_solution.view(), + other_solution = other_solution.view(), + delta_vector = make_span(delta_vector), + n_points_to_search] __device__(i_t idx) { + f_t guiding_val = guiding_solution.assignment[idx]; + f_t other_val = other_solution.assignment[idx]; + delta_vector[idx] = (other_val - guiding_val) / (n_points_to_search + 1); + }); + return delta_vector; + } + std::pair, bool> recombine(solution_t& a, solution_t& b, const weight_t& weights) { raft::common::nvtx::range fun_scope("line_segment_recombiner"); + auto& guiding_solution = a.get_feasible() ? a : b; + auto& other_solution = a.get_feasible() ? b : a; // copy the solution from A - solution_t offspring(a); - // TODO test the time limit of two seconds - timer_t line_segment_timer{2.}; + solution_t offspring(guiding_solution); + timer_t line_segment_timer{ls_recombiner_config_t::time_limit}; // TODO after we have the conic combination, detect the lambda change // (i.e. the integral variables flip on line segment) - i_t n_points_to_search = 20; - bool is_feasibility_run = false; + i_t n_points_to_search = ls_recombiner_config_t::n_points_to_search; + const bool is_feasibility_run = false; + i_t n_different_vars = + this->assign_same_integer_values(guiding_solution, other_solution, offspring); + rmm::device_uvector delta_vector = generate_delta_vector( + guiding_solution, other_solution, offspring, n_points_to_search, n_different_vars); line_segment_search.fj.copy_weights(weights, offspring.handle_ptr); + // return infeasible results as well + line_segment_search.settings.recombiner_mode = true; + line_segment_search.settings.best_of_parents_cost = + std::min(guiding_solution.get_quality(weights), other_solution.get_quality(weights)); + line_segment_search.settings.parents_infeasible = + !guiding_solution.get_feasible() && !other_solution.get_feasible(); // TODO fix common part and run FJ on remaining line_segment_search.search_line_segment(offspring, - a.assignment, - b.assignment, + guiding_solution.assignment, + other_solution.assignment, n_points_to_search, + delta_vector, is_feasibility_run, line_segment_timer); - bool same_as_parents = this->check_if_offspring_is_same_as_parents(offspring, a, b); + line_segment_search.settings = {}; + bool better_cost_than_parents = + offspring.get_quality(weights) < + std::min(other_solution.get_quality(weights), guiding_solution.get_quality(weights)); + bool better_feasibility_than_parents = offspring.get_feasible() && + !other_solution.get_feasible() && + !guiding_solution.get_feasible(); + bool same_as_parents = + this->check_if_offspring_is_same_as_parents(offspring, guiding_solution, other_solution); + // adjust the max_n_of_vars_from_other + if (n_different_vars > (i_t)ls_recombiner_config_t::max_n_of_vars_from_other) { + if (same_as_parents) { + ls_recombiner_config_t::increase_max_n_of_vars_from_other(); + } else { + ls_recombiner_config_t::decrease_max_n_of_vars_from_other(); + } + } + if (better_cost_than_parents || better_feasibility_than_parents) { + CUOPT_LOG_DEBUG("Offspring is feasible or better than both parents"); + return std::make_pair(offspring, true); + } return std::make_pair(offspring, !same_as_parents); } diff --git a/cpp/src/mip/diversity/recombiners/recombiner.cuh b/cpp/src/mip/diversity/recombiners/recombiner.cuh index 2f1243f1f..e4a86b491 100644 --- a/cpp/src/mip/diversity/recombiners/recombiner.cuh +++ b/cpp/src/mip/diversity/recombiners/recombiner.cuh @@ -17,14 +17,19 @@ #pragma once +#include "recombiner_configs.hpp" + #include #include #include #include #include +#include #include +#include #include +#include namespace cuopt::linear_programming::detail { @@ -50,16 +55,11 @@ __global__ void assign_same_variables_kernel(typename solution_t::view raft::device_span remaining_indices, i_t* n_remaining) { - if (TH_ID_X >= a.assignment.size()) return; - const i_t var_idx = TH_ID_X; - bool is_integer_var = a.problem.is_integer_var(var_idx); - - if (diverse_equal(a.assignment[var_idx], - b.assignment[var_idx], - a.problem.variable_lower_bounds[var_idx], - a.problem.variable_upper_bounds[var_idx], - is_integer_var, - a.problem.tolerances.integrality_tolerance)) { + if (TH_ID_X >= a.problem.n_integer_vars) return; + const i_t var_idx = a.problem.integer_indices[TH_ID_X]; + + if (integer_equal( + a.assignment[var_idx], b.assignment[var_idx], a.problem.tolerances.integrality_tolerance)) { offspring.assignment[var_idx] = a.assignment[var_idx]; } else { i_t idx = atomicAdd(n_remaining, 1); @@ -71,33 +71,79 @@ template class recombiner_t { public: recombiner_t(mip_solver_context_t& context_, - i_t n_vars, + i_t n_integer_vars, const raft::handle_t* handle_ptr) : context(context_), - remaining_indices(n_vars, handle_ptr->get_stream()), + remaining_indices(n_integer_vars, handle_ptr->get_stream()), n_remaining(handle_ptr->get_stream()) { } - void reset(i_t n_vars, const raft::handle_t* handle_ptr) + void reset(i_t n_integer_vars, const raft::handle_t* handle_ptr) { n_remaining.set_value_to_zero_async(handle_ptr->get_stream()); - remaining_indices.resize(n_vars, handle_ptr->get_stream()); + remaining_indices.resize(n_integer_vars, handle_ptr->get_stream()); } - void assign_same_integer_values(solution_t& a, - solution_t& b, - solution_t& offspring) + i_t assign_same_integer_values(solution_t& a, + solution_t& b, + solution_t& offspring) { - reset(a.problem_ptr->n_variables, a.handle_ptr); + reset(a.problem_ptr->n_integer_vars, a.handle_ptr); const i_t TPB = 128; - i_t n_blocks = (a.problem_ptr->n_variables + TPB - 1) / TPB; + i_t n_blocks = (a.problem_ptr->n_integer_vars + TPB - 1) / TPB; assign_same_variables_kernel <<get_stream()>>>(a.view(), b.view(), offspring.view(), cuopt::make_span(remaining_indices), n_remaining.data()); + i_t remaining_variables = this->n_remaining.value(a.handle_ptr->get_stream()); + + auto vec_remaining_indices = + host_copy(this->remaining_indices.data(), remaining_variables, a.handle_ptr->get_stream()); + auto vec_objective_coeffs = host_copy(offspring.problem_ptr->objective_coefficients.data(), + offspring.problem_ptr->n_variables, + a.handle_ptr->get_stream()); + auto integer_indices = host_copy(offspring.problem_ptr->integer_indices.data(), + offspring.problem_ptr->n_integer_vars, + a.handle_ptr->get_stream()); + std::vector objective_indices; + for (size_t i = 0; i < vec_objective_coeffs.size(); ++i) { + if (vec_objective_coeffs[i] != 0 && + std::find(integer_indices.begin(), integer_indices.end(), i) != integer_indices.end()) { + objective_indices.push_back(i); + } + } + std::vector objective_indices_in_subproblem; + for (auto var : vec_remaining_indices) { + if (vec_objective_coeffs[var] != 0) { objective_indices_in_subproblem.push_back(var); } + } + CUOPT_LOG_DEBUG("n_objective_vars in different vars %d n_objective_vars %d", + objective_indices_in_subproblem.size(), + objective_indices.size()); + if (objective_indices_in_subproblem.size() < 0.4 * remaining_variables) { + std::default_random_engine rng_host(cuopt::seed_generator::get_seed()); + std::vector objective_indices_not_in_subproblem; + std::set_difference(objective_indices.begin(), + objective_indices.end(), + objective_indices_in_subproblem.begin(), + objective_indices_in_subproblem.end(), + std::back_inserter(objective_indices_not_in_subproblem)); + std::shuffle(objective_indices_not_in_subproblem.begin(), + objective_indices_not_in_subproblem.end(), + rng_host); + for (auto var : objective_indices_not_in_subproblem) { + objective_indices_in_subproblem.push_back(var); + vec_remaining_indices.push_back(var); + if (objective_indices_in_subproblem.size() >= 0.4 * remaining_variables) { break; } + } + } + raft::copy(this->remaining_indices.data(), + vec_remaining_indices.data(), + vec_remaining_indices.size(), + a.handle_ptr->get_stream()); + return vec_remaining_indices.size(); } bool check_if_offspring_is_same_as_parents(solution_t& offspring, @@ -109,13 +155,48 @@ class recombiner_t { a.assignment, a.problem_ptr->tolerances.integrality_tolerance, a.handle_ptr); - if (equals_a) { return true; } + if (equals_a) { + CUOPT_LOG_DEBUG("Offspring is same as parent guiding!"); + return true; + } bool equals_b = check_integer_equal_on_indices(offspring.problem_ptr->integer_indices, offspring.assignment, b.assignment, b.problem_ptr->tolerances.integrality_tolerance, b.handle_ptr); - return equals_b; + if (equals_b) { + CUOPT_LOG_DEBUG("Offspring is same as parent other!"); + return true; + } + return false; + } + + void compute_vars_to_fix(solution_t& offspring, + rmm::device_uvector& vars_to_fix, + i_t n_vars_from_other, + i_t n_vars_from_guiding) + { + vars_to_fix.resize(n_vars_from_guiding, offspring.handle_ptr->get_stream()); + // set difference needs two sorted arrays + thrust::sort(offspring.handle_ptr->get_thrust_policy(), + this->remaining_indices.data(), + this->remaining_indices.data() + n_vars_from_other); + cuopt_assert((thrust::is_sorted(offspring.handle_ptr->get_thrust_policy(), + offspring.problem_ptr->integer_indices.begin(), + offspring.problem_ptr->integer_indices.end())), + "vars_to_fix should be sorted!"); + // get the variables to fix (common variables) + auto iter = thrust::set_difference(offspring.handle_ptr->get_thrust_policy(), + offspring.problem_ptr->integer_indices.begin(), + offspring.problem_ptr->integer_indices.end(), + this->remaining_indices.data(), + this->remaining_indices.data() + n_vars_from_other, + vars_to_fix.begin()); + cuopt_assert(iter - vars_to_fix.begin() == n_vars_from_guiding, "The size should match!"); + cuopt_assert((thrust::is_sorted(offspring.handle_ptr->get_thrust_policy(), + vars_to_fix.data(), + vars_to_fix.data() + n_vars_from_guiding)), + "vars_to_fix should be sorted!"); } mip_solver_context_t& context; diff --git a/cpp/src/mip/diversity/recombiners/recombiner_configs.hpp b/cpp/src/mip/diversity/recombiners/recombiner_configs.hpp new file mode 100644 index 000000000..c4c3ed0b9 --- /dev/null +++ b/cpp/src/mip/diversity/recombiners/recombiner_configs.hpp @@ -0,0 +1,104 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +namespace cuopt::linear_programming::detail { + +struct bp_recombiner_config_t { + static constexpr double bounds_prop_time_limit = 2; + static constexpr double lp_after_bounds_prop_time_limit = 2; + // number of repair iterations even if it fails during the repair + static constexpr size_t n_repair_iterations = 10; + static constexpr size_t initial_n_of_vars_from_other = 200; + static constexpr size_t max_different_var_limit = 10000; + static constexpr size_t min_different_var_limit = 20; + static size_t max_n_of_vars_from_other; + static constexpr double n_var_ratio_increase_factor = 1.1; + static constexpr double n_var_ratio_decrease_factor = 0.99; + static void increase_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::min( + size_t(max_n_of_vars_from_other * n_var_ratio_increase_factor), max_different_var_limit); + CUOPT_LOG_DEBUG("Increased max_n_of_vars_from_other in BP recombiner to %lu", + max_n_of_vars_from_other); + } + static void decrease_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::max( + size_t(max_n_of_vars_from_other * n_var_ratio_decrease_factor), min_different_var_limit); + CUOPT_LOG_DEBUG("Decreased max_n_of_vars_from_other in BP recombiner to %lu", + max_n_of_vars_from_other); + } +}; + +struct ls_recombiner_config_t { + // line segment related configs + // FIXME: not implemented yet + static constexpr bool use_fj_for_rounding = false; + static constexpr int n_points_to_search = 20; + static constexpr double time_limit = 2.; + static constexpr size_t initial_n_of_vars_from_other = 200; + static constexpr size_t max_different_var_limit = 10000; + static constexpr size_t min_different_var_limit = 20; + static size_t max_n_of_vars_from_other; + static constexpr double n_var_ratio_increase_factor = 1.1; + static constexpr double n_var_ratio_decrease_factor = 0.99; + static void increase_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::min( + size_t(max_n_of_vars_from_other * n_var_ratio_increase_factor), max_different_var_limit); + CUOPT_LOG_DEBUG("Increased max_n_of_vars_from_other in LS recombiner to %lu", + max_n_of_vars_from_other); + } + static void decrease_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::max( + size_t(max_n_of_vars_from_other * n_var_ratio_decrease_factor), min_different_var_limit); + CUOPT_LOG_DEBUG("Decreased max_n_of_vars_from_other in LS recombiner to %lu", + max_n_of_vars_from_other); + } +}; + +struct fp_recombiner_config_t { + static constexpr double infeasibility_detection_time_limit = 0.5; + static constexpr double fp_time_limit = 3.; + static constexpr double alpha = 0.99; + static constexpr double alpha_decrease_factor = 0.9; + static constexpr size_t initial_n_of_vars_from_other = 200; + static constexpr size_t max_different_var_limit = 10000; + static constexpr size_t min_different_var_limit = 20; + static size_t max_n_of_vars_from_other; + static constexpr double n_var_ratio_increase_factor = 1.1; + static constexpr double n_var_ratio_decrease_factor = 0.99; + static void increase_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::min( + size_t(max_n_of_vars_from_other * n_var_ratio_increase_factor), max_different_var_limit); + CUOPT_LOG_DEBUG("Increased max_n_of_vars_from_other in FP recombiner to %lu", + max_n_of_vars_from_other); + } + static void decrease_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::max( + size_t(max_n_of_vars_from_other * n_var_ratio_decrease_factor), min_different_var_limit); + CUOPT_LOG_DEBUG("Decreased max_n_of_vars_from_other in FP recombiner to %lu", + max_n_of_vars_from_other); + } +}; + +} // namespace cuopt::linear_programming::detail \ No newline at end of file diff --git a/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp b/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp index 36091db66..1c34b0b5c 100644 --- a/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp +++ b/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp @@ -37,10 +37,15 @@ struct recombine_stats { void add_success() { ++success; } - void update_improve_stats(double cost_new, double cost_first, double cost_second) + bool update_improve_stats(double cost_new, double cost_first, double cost_second) { - if (cost_new < (std::min(cost_first, cost_second) - OBJECTIVE_EPSILON)) ++better_than_both; + bool is_better_than_both = false; + if (cost_new < (std::min(cost_first, cost_second) - OBJECTIVE_EPSILON)) { + ++better_than_both; + is_better_than_both = true; + } if (cost_new < (std::max(cost_first, cost_second) - OBJECTIVE_EPSILON)) ++better_than_one; + return is_better_than_both; } void add_attempt() { ++attempts; } @@ -67,6 +72,21 @@ struct all_recombine_stats { // enum of the last attempted recombiner std::optional last_attempt; + double last_recombiner_time; + std::chrono::high_resolution_clock::time_point last_recombiner_start_time; + + void start_recombiner_time() + { + last_recombiner_start_time = std::chrono::high_resolution_clock::now(); + } + void stop_recombiner_time() + { + last_recombiner_time = std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - last_recombiner_start_time) + .count(); + } + + double get_last_recombiner_time() { return last_recombiner_time; } void reset() { @@ -76,6 +96,8 @@ struct all_recombine_stats { last_attempt.reset(); } + recombiner_enum_t get_last_attempt() { return last_attempt.value(); } + void add_attempt(recombiner_enum_t r) { last_attempt = r; @@ -84,9 +106,9 @@ struct all_recombine_stats { void add_success() { stats[static_cast(last_attempt.value())].add_success(); } - void update_improve_stats(double cost_new, double cost_first, double cost_second) + bool update_improve_stats(double cost_new, double cost_first, double cost_second) { - stats[static_cast(last_attempt.value())].update_improve_stats( + return stats[static_cast(last_attempt.value())].update_improve_stats( cost_new, cost_first, cost_second); } diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump.cu b/cpp/src/mip/feasibility_jump/feasibility_jump.cu index 724a053cf..5043d7780 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump.cu +++ b/cpp/src/mip/feasibility_jump/feasibility_jump.cu @@ -30,6 +30,7 @@ #include #include +#include #include #include #include @@ -158,7 +159,6 @@ fj_t::climber_data_t::view_t fj_t::climber_data_t::view() v.iteration_related_variables = iteration_related_variables.view(); v.constraints_changed = make_span(constraints_changed); v.best_assignment = make_span(best_assignment); - v.farthest_l1_sol = make_span(farthest_l1_sol); v.incumbent_assignment = make_span(incumbent_assignment); v.cstr_weights = make_span(fj.cstr_weights); v.cstr_left_weights = make_span(fj.cstr_left_weights); @@ -196,6 +196,9 @@ fj_t::climber_data_t::view_t fj_t::climber_data_t::view() v.work_id_to_bin_var_idx = make_span(fj.work_id_to_bin_var_idx); v.work_id_to_nonbin_var_idx = make_span(fj.work_id_to_nonbin_var_idx); v.work_ids_for_related_vars = make_span(fj.work_ids_for_related_vars); + v.fractional_variables = fractional_variables.view(); + v.saved_best_fractional_count = saved_best_fractional_count.data(); + v.handle_fractionals_only = handle_fractionals_only.data(); v.selected_var = selected_var.data(); v.violation_score = violation_score.data(); v.weighted_violation_score = weighted_violation_score.data(); @@ -203,7 +206,6 @@ fj_t::climber_data_t::view_t fj_t::climber_data_t::view() v.local_minimums_reached = local_minimums_reached.data(); v.iterations = iterations.data(); v.best_excess = best_excess.data(); - v.best_l1_distance = best_l1_distance.data(); v.best_objective = best_objective.data(); v.saved_solution_objective = saved_solution_objective.data(); v.incumbent_quality = incumbent_quality.data(); @@ -218,6 +220,7 @@ fj_t::climber_data_t::view_t fj_t::climber_data_t::view() v.break_condition = break_condition.data(); v.temp_break_condition = temp_break_condition.data(); v.best_jump_idx = best_jump_idx.data(); + v.small_move_tabu = small_move_tabu.data(); v.stop_threshold = fj.stop_threshold; v.iterations_until_feasible_counter = iterations_until_feasible_counter.data(); v.full_refresh_iteration = full_refresh_iteration.data(); @@ -256,6 +259,7 @@ void fj_t::climber_data_t::clear_sets(const rmm::cuda_stream_view& str violated_constraints.clear(stream); candidate_variables.clear(stream); iteration_related_variables.clear(stream); + fractional_variables.clear(stream); } template @@ -363,9 +367,10 @@ void fj_t::climber_init(i_t climber_idx, const rmm::cuda_stream_view& thrust::counting_iterator(0), thrust::counting_iterator(pb_ptr->n_variables), [pb = pb_ptr->view(), + mode = settings.mode, incumbent_assignment = climber->incumbent_assignment.data()] __device__(i_t var_idx) { // round if integer - if (pb.is_integer_var(var_idx)) { + if (mode != fj_mode_t::ROUNDING && pb.is_integer_var(var_idx)) { incumbent_assignment[var_idx] = round(incumbent_assignment[var_idx]); } // clamp to bounds @@ -374,6 +379,20 @@ void fj_t::climber_init(i_t climber_idx, const rmm::cuda_stream_view& min(pb.variable_upper_bounds[var_idx], incumbent_assignment[var_idx])); }); + thrust::for_each( + rmm::exec_policy(climber_stream), + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(pb_ptr->n_variables), + [v = view] __device__(i_t var_idx) { + if (v.pb.is_integer_var(var_idx) && !v.pb.is_integer(v.incumbent_assignment[var_idx])) + v.fractional_variables.insert(var_idx); + }); + + i_t fractional_var_count = climber->fractional_variables.set_size.value(climber_stream); + climber->saved_best_fractional_count.set_value_async(fractional_var_count, climber_stream); + climber->handle_fractionals_only.set_value_to_zero_async(climber_stream); + CUOPT_LOG_TRACE("fractional_var_count = %d\n", fractional_var_count); + objective_vars.resize(pb_ptr->n_variables, climber_stream); auto end = thrust::copy_if(rmm::exec_policy(climber_stream), thrust::counting_iterator(0), @@ -392,14 +411,14 @@ void fj_t::climber_init(i_t climber_idx, const rmm::cuda_stream_view& climber->best_objective.set_value_async(inf, climber_stream); climber->saved_solution_objective.set_value_async(inf, climber_stream); climber->violation_score.set_value_to_zero_async(climber_stream); - climber->best_l1_distance.set_value_to_zero_async(climber_stream); climber->weighted_violation_score.set_value_to_zero_async(climber_stream); init_lhs_and_violation<<<256, 256, 0, climber_stream.value()>>>(view); // initialize the best_objective values according to the initial assignment f_t best_obj = compute_objective_from_vec( climber->incumbent_assignment, pb_ptr->objective_coefficients, climber_stream); - if (climber->violated_constraints.set_size.value(climber_stream) == 0) { + if (climber->violated_constraints.set_size.value(climber_stream) == 0 && + (settings.mode != fj_mode_t::ROUNDING || fractional_var_count == 0)) { climber->best_excess.set_value_to_zero_async(climber_stream); climber->best_objective.set_value_async(best_obj, climber_stream); climber->saved_solution_objective.set_value_async(best_obj, climber_stream); @@ -421,6 +440,9 @@ void fj_t::climber_init(i_t climber_idx, const rmm::cuda_stream_view& climber->iterations.set_value_to_zero_async(climber_stream); climber->full_refresh_iteration.set_value_to_zero_async(climber_stream); climber->iterations_until_feasible_counter.set_value_to_zero_async(climber_stream); + climber->small_move_tabu.set_value_to_zero_async(climber_stream); + + climber_stream.synchronize(); climber_stream.synchronize(); @@ -615,10 +637,14 @@ void fj_t::run_step_device(const rmm::cuda_stream_view& climber_stream // ensure an updated copy of the settings is used device-side raft::copy(v.settings, &settings, 1, climber_stream); - bool is_binary_pb = pb_ptr->n_variables == thrust::count(handle_ptr->get_thrust_policy(), + bool is_binary_pb = pb_ptr->n_variables == thrust::count(handle_ptr->get_thrust_policy(), pb_ptr->is_binary_variable.begin(), pb_ptr->is_binary_variable.end(), 1); + // if we're in rounding mode, do not treat the problem as a purely binary one + // as it breaks assumptions in the binary_pb codepath + if (settings.mode == fj_mode_t::ROUNDING) { is_binary_pb = false; } + bool use_load_balancing = false; if (settings.load_balancing_mode == fj_load_balancing_mode_t::ALWAYS_OFF) { use_load_balancing = false; @@ -628,6 +654,8 @@ void fj_t::run_step_device(const rmm::cuda_stream_view& climber_stream use_load_balancing = pb_ptr->n_variables > settings.parameters.load_balancing_codepath_min_varcount; } + // Load-balanced codepath not updated yet to handle rounding mode + if (settings.mode == fj_mode_t::ROUNDING) { use_load_balancing = false; } cudaGraph_t graph; void* kernel_args[] = {&v}; @@ -767,6 +795,24 @@ void fj_t::run_step_device(const rmm::cuda_stream_view& climber_stream if (use_graph) cudaGraphLaunch(graph_instance, climber_stream); } +template +void fj_t::round_remaining_fractionals(solution_t& solution, i_t climber_idx) +{ + auto& data = *climbers[climber_idx]; + + auto climber_stream = data.stream.view(); + if (climber_idx == 0) climber_stream = handle_ptr->get_stream(); + + bool handle_fractionals_only = true; + data.handle_fractionals_only.set_value_async(handle_fractionals_only, climber_stream); + data.break_condition.set_value_to_zero_async(climber_stream); + data.temp_break_condition.set_value_to_zero_async(climber_stream); + climber_stream.synchronize(); + + // Run the fractional move selection and assignment update kernels until all have been rounded + host_loop(solution, climber_idx); +} + template void fj_t::refresh_lhs_and_violation(const rmm::cuda_stream_view& stream, i_t climber_idx) { @@ -894,11 +940,15 @@ i_t fj_t::host_loop(solution_t& solution, i_t climber_idx) same_sol, max_cstr_weight.value(climber_stream)); #endif - CUOPT_LOG_TRACE("EXIT step %d, best objective %f best_excess %f, feas %d", + CUOPT_LOG_DEBUG("EXIT FJ step %d, best objective %f best_excess %f, feas %d, local mins %d", data.iterations.value(climber_stream), solution.get_user_objective(), solution.get_total_excess(), - solution.get_feasible()); + solution.get_feasible(), + data.local_minimums_reached.value(climber_stream)); + + CUOPT_LOG_TRACE("best fractional count %d", + data.saved_best_fractional_count.value(climber_stream)); return steps; } @@ -942,7 +992,6 @@ void fj_t::resize_vectors(const raft::handle_t* handle_ptr) climbers[0]->constraints_changed.resize(pb_ptr->n_constraints, handle_ptr->get_stream()); climbers[0]->violated_constraints.resize(pb_ptr->n_constraints, handle_ptr->get_stream()); climbers[0]->best_assignment.resize(pb_ptr->n_variables, handle_ptr->get_stream()); - climbers[0]->farthest_l1_sol.resize(pb_ptr->n_variables, handle_ptr->get_stream()); climbers[0]->incumbent_assignment.resize(pb_ptr->n_variables, handle_ptr->get_stream()); climbers[0]->incumbent_lhs.resize(pb_ptr->n_constraints, handle_ptr->get_stream()); climbers[0]->incumbent_lhs_sumcomp.resize(pb_ptr->n_constraints, handle_ptr->get_stream()); @@ -990,12 +1039,24 @@ i_t fj_t::solve(solution_t& solution) timer_t timer(settings.time_limit); handle_ptr = const_cast(solution.handle_ptr); pb_ptr = solution.problem_ptr; - cuopt_func_call(solution.test_variable_bounds(true)); - cuopt_assert(solution.test_number_all_integer(), "All integers must be rounded"); + if (settings.mode != fj_mode_t::ROUNDING) { + cuopt_func_call(solution.test_variable_bounds(true)); + cuopt_assert(solution.test_number_all_integer(), "All integers must be rounded"); + } pb_ptr->check_problem_representation(true); resize_vectors(solution.handle_ptr); bool is_initial_feasible = solution.compute_feasibility(); + // if we're in rounding mode, split the time/iteration limit between the first and second stage + cuopt_assert(settings.parameters.rounding_second_stage_split >= 0 && + settings.parameters.rounding_second_stage_split <= 1, + "rounding_second_stage_split must be between 0 and 1"); + if (settings.mode == fj_mode_t::ROUNDING) { + settings.time_limit = + settings.time_limit * (1 - settings.parameters.rounding_second_stage_split); + settings.iteration_limit = + settings.iteration_limit * (1 - settings.parameters.rounding_second_stage_split); + } // TODO only call this when the size is different device_init(handle_ptr->get_stream()); @@ -1024,9 +1085,35 @@ i_t fj_t::solve(solution_t& solution) f_t effort_rate = (f_t)iterations / timer.elapsed_time(); + // If we're in rounding mode and some fractionals remain: round them all + // limit = total_limit * second_stage_split + if (settings.mode == fj_mode_t::ROUNDING && + climbers[0]->fractional_variables.set_size.value(handle_ptr->get_stream()) > 0) { + settings.time_limit = settings.time_limit * settings.parameters.rounding_second_stage_split; + settings.iteration_limit = + settings.iteration_limit * settings.parameters.rounding_second_stage_split; + + round_remaining_fractionals(solution); + // if time limit exceeded: round all remaining fractionals if any by nearest rounding. + if (climbers[0]->fractional_variables.set_size.value(handle_ptr->get_stream()) > 0) { + solution.round_nearest(); + } + } + CUOPT_LOG_TRACE("GPU solver took %g", timer.elapsed_time()); CUOPT_LOG_TRACE("limit reached, effort rate %g steps/secm %d steps", effort_rate, iterations); reset_cuda_graph(); + i_t n_integer_vars = thrust::count_if( + handle_ptr->get_thrust_policy(), + solution.problem_ptr->integer_indices.begin(), + solution.problem_ptr->integer_indices.end(), + [pb = solution.problem_ptr->view(), + assignment_ptr = solution.assignment.data()] __device__(i_t idx) -> bool { + if (!pb.is_integer(assignment_ptr[idx])) { + DEVICE_LOG_ERROR("variable %d is not integer, value %g\n", idx, assignment_ptr[idx]); + } + return pb.is_integer(assignment_ptr[idx]); + }); cuopt_assert(solution.test_number_all_integer(), "All integers must be rounded"); bool is_new_feasible = solution.compute_feasibility(); diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump.cuh b/cpp/src/mip/feasibility_jump/feasibility_jump.cuh index 673229d4d..a55f115f1 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump.cuh +++ b/cpp/src/mip/feasibility_jump/feasibility_jump.cuh @@ -71,6 +71,12 @@ struct fj_hyper_parameters_t { double excess_improvement_weight = (1.0 / 2.0); double weight_smoothing_probability = 0.0003; + double fractional_score_multiplier = 100; + double rounding_second_stage_split = 0.1; + + double small_move_tabu_threshold = 1e-6; + int small_move_tabu_tenure = 4; + // load-balancing related settings int old_codepath_total_var_to_relvar_ratio_threshold = 200; int load_balancing_codepath_min_varcount = 3200; @@ -94,6 +100,7 @@ enum class fj_mode_t { FIRST_FEASIBLE, // iterate until a feasible solution is found, then return GREEDY_DESCENT, // single descent until no improving jumps can be made TREE, // tree mode + ROUNDING, // FJ as rounding procedure for fractionals EXIT_NON_IMPROVING // iterate until we are don't improve the best }; @@ -113,8 +120,8 @@ struct fj_settings_t { double infeasibility_weight = 1.0; bool update_weights = true; bool feasibility_run = true; - bool keep_farthest_l1_sol = false; fj_load_balancing_mode_t load_balancing_mode{fj_load_balancing_mode_t::AUTO}; + double baseline_objective_for_longer_run{std::numeric_limits::lowest()}; }; struct fj_move_t { @@ -215,6 +222,8 @@ class fj_t { void refresh_lhs_and_violation(const rmm::cuda_stream_view& stream, i_t climber_idx = 0); // load balancing void load_balancing_score_update(const rmm::cuda_stream_view& stream, i_t climber_idx = 0); + // executed after a roudning FJ run if any fractionals remain to eliminate them + void round_remaining_fractionals(solution_t& solution, i_t climber_idx = 0); public: mip_solver_context_t& context; @@ -284,7 +293,6 @@ class fj_t { rmm::device_scalar local_minimums_reached; rmm::device_scalar iterations; rmm::device_scalar best_excess; - rmm::device_scalar best_l1_distance; rmm::device_scalar best_objective; rmm::device_scalar saved_solution_objective; rmm::device_scalar incumbent_quality; @@ -308,7 +316,6 @@ class fj_t { bitmap_t iteration_related_variables; rmm::device_uvector constraints_changed; rmm::device_uvector best_assignment; - rmm::device_uvector farthest_l1_sol; rmm::device_uvector incumbent_assignment; rmm::device_uvector incumbent_lhs; // compensation term of the Kahan summation algorithm @@ -328,6 +335,12 @@ class fj_t { rmm::device_uvector jump_candidates; rmm::device_uvector jump_candidate_count; rmm::device_uvector jump_locks; + rmm::device_scalar small_move_tabu; + + // ROUNDING mode related members + contiguous_set_t fractional_variables; + rmm::device_scalar saved_best_fractional_count; + rmm::device_scalar handle_fractionals_only; rmm::device_uvector grid_score_buf; rmm::device_uvector grid_var_buf; @@ -347,7 +360,6 @@ class fj_t { local_minimums_reached(0, fj.handle_ptr->get_stream()), iterations(0, fj.handle_ptr->get_stream()), best_excess(-std::numeric_limits::infinity(), fj.handle_ptr->get_stream()), - best_l1_distance(0., fj.handle_ptr->get_stream()), best_objective(+std::numeric_limits::infinity(), fj.handle_ptr->get_stream()), saved_solution_objective(+std::numeric_limits::infinity(), fj.handle_ptr->get_stream()), @@ -361,7 +373,6 @@ class fj_t { iteration_related_variables(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), constraints_changed(fj.pb_ptr->n_constraints, fj.handle_ptr->get_stream()), best_assignment(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), - farthest_l1_sol(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), tabu_nodec_until(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), tabu_noinc_until(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), tabu_lastdec(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), @@ -380,6 +391,10 @@ class fj_t { jump_candidates(fj.pb_ptr->coefficients.size(), fj.handle_ptr->get_stream()), jump_candidate_count(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), jump_locks(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), + fractional_variables(fj.pb_ptr->n_variables, fj.handle_ptr->get_stream()), + small_move_tabu(0, fj.handle_ptr->get_stream()), + handle_fractionals_only(false, fj.handle_ptr->get_stream()), + saved_best_fractional_count(0, fj.handle_ptr->get_stream()), candidate_arrived_workids(fj.pb_ptr->coefficients.size(), fj.handle_ptr->get_stream()), grid_score_buf(0, fj.handle_ptr->get_stream()), grid_var_buf(0, fj.handle_ptr->get_stream()), @@ -416,7 +431,6 @@ class fj_t { raft::device_span cstr_left_weights; raft::device_span incumbent_assignment; raft::device_span best_assignment; - raft::device_span farthest_l1_sol; raft::device_span incumbent_lhs; raft::device_span incumbent_lhs_sumcomp; raft::device_span jump_move_scores; @@ -451,6 +465,10 @@ class fj_t { typename contiguous_set_t::view_t candidate_variables; typename bitmap_t::view_t iteration_related_variables; + // ROUNDING mode related members + typename contiguous_set_t::view_t fractional_variables; + i_t* saved_best_fractional_count; + bool* handle_fractionals_only; // load balancing structures raft::device_span row_size_bin_prefix_sum; raft::device_span row_size_nonbin_prefix_sum; @@ -465,7 +483,6 @@ class fj_t { i_t* iterations; i_t* last_iter_candidates; f_t* best_excess; - f_t* best_l1_distance; f_t* best_objective; f_t* saved_solution_objective; f_t* incumbent_quality; @@ -476,6 +493,7 @@ class fj_t { i_t* full_refresh_iteration; cub::KeyValuePair* best_jump_idx; f_t stop_threshold; + i_t* small_move_tabu; i_t* last_minimum_iteration; i_t* last_improving_minimum; @@ -546,6 +564,8 @@ class fj_t { (delta > 0 && iter < tabu_noinc_until[var_idx])) return false; + if (*handle_fractionals_only && !fractional_variables.contains(var_idx)) return false; + // give priority to MTM moves. if (jump_move_scores[var_idx].base > 0) return true; diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu index 7ff901ca7..213d31bc5 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu +++ b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu @@ -115,6 +115,7 @@ DI std::pair feas_score_constraint( f_t new_slack = rhs - new_lhs; cuopt_assert(isfinite(cstr_weight), "invalid weight"); + cuopt_assert(cstr_weight >= 0, "invalid weight"); cuopt_assert(isfinite(old_lhs), ""); cuopt_assert(isfinite(new_lhs), ""); cuopt_assert(isfinite(old_slack) && isfinite(new_slack), ""); @@ -339,7 +340,16 @@ DI typename fj_t::move_score_info_t compute_new_score( typename fj_t::move_score_info_t score_info; - score_info.score.base = round(base_obj + base_feas_sum); + // Add a large bonus if this move would turn a fractional variable integral + // in rounding mode + f_t integrality_multiplier = 1; + if (fj.settings->mode == fj_mode_t::ROUNDING) { + if (fj.pb.is_integer_var(var_idx) && !fj.pb.is_integer(fj.incumbent_assignment[var_idx]) && + fj.pb.is_integer(fj.incumbent_assignment[var_idx] + delta)) { + integrality_multiplier = fj.settings->parameters.fractional_score_multiplier; + } + } + score_info.score.base = round(base_obj + base_feas_sum) * integrality_multiplier; score_info.score.bonus = round(bonus_breakthrough + bonus_robust_sum); score_info.infeasibility = base_feas_sum; @@ -352,6 +362,29 @@ DI typename fj_t::move_score_info_t compute_new_score( return score_info; } +// Returns the current slack, and the variable delta that would nullify this slack ("tighten" it) +template +DI thrust::tuple get_mtm_for_bound( + const typename fj_t::climber_data_t::view_t& fj, + i_t var_idx, + i_t cstr_idx, + f_t cstr_coeff, + f_t bound, + f_t sign) +{ + f_t delta_ij = 0; + f_t slack = 0; + f_t old_val = fj.incumbent_assignment[var_idx]; + + f_t lhs = fj.incumbent_lhs[cstr_idx] * sign; + f_t rhs = bound * sign; + slack = rhs - lhs; // bound might be infinite. let the caller handle this case + + delta_ij = slack / (cstr_coeff * sign); + + return {delta_ij, slack}; +} + template DI thrust::tuple get_mtm_for_constraint( const typename fj_t::climber_data_t::view_t& fj, @@ -435,6 +468,9 @@ DI std::pair::move_score_info_t> compute_best_mtm( // reservoir sampling to get 20 random constraints f_t selection_probability = min(1.0, 20.0 / (offset_end - offset_begin)); if (rng.next_double() > selection_probability) continue; + } else if constexpr (move_type == MTMMoveType::FJ_MTM_ALL) { + // sample all constraints, regardless of violation status + ; } f_t new_val; @@ -457,6 +493,8 @@ DI std::pair::move_score_info_t> compute_best_mtm( } } + if (fj.pb.is_integer_var(var_idx)) new_val = round(new_val); + if (fj.pb.integer_equal(new_val, old_val) || !isfinite(new_val)) continue; if (threadIdx.x == 0) { @@ -506,9 +544,10 @@ DI void update_jump_value(typename fj_t::climber_data_t::view_t fj, i_ auto best_score_info = fj_t::move_score_info_t::invalid(); if constexpr (!is_binary_pb) { if (fj.pb.is_binary_variable[var_idx] && fj.pb.is_integer(fj.incumbent_assignment[var_idx])) { - delta = 1.0 - 2 * fj.incumbent_assignment[var_idx]; + delta = round(1.0 - 2 * fj.incumbent_assignment[var_idx]); if (threadIdx.x == 0) { - cuopt_assert(fj.incumbent_assignment[var_idx] == 0 || fj.incumbent_assignment[var_idx] == 1, + cuopt_assert(fj.pb.integer_equal(fj.incumbent_assignment[var_idx], 0) || + fj.pb.integer_equal(fj.incumbent_assignment[var_idx], 1), "Current assignment is not binary!"); cuopt_assert( fj.pb.variable_lower_bounds[var_idx] == 0 && fj.pb.variable_upper_bounds[var_idx] == 1, @@ -525,7 +564,7 @@ DI void update_jump_value(typename fj_t::climber_data_t::view_t fj, i_ best_score_info = score_info; } } else { - delta = 1.0 - 2 * fj.incumbent_assignment[var_idx]; + delta = round(1.0 - 2 * fj.incumbent_assignment[var_idx]); if (threadIdx.x == 0) { cuopt_assert( fj.pb.variable_lower_bounds[var_idx] == 0 && fj.pb.variable_upper_bounds[var_idx] == 1, ""); @@ -579,34 +618,23 @@ DI bool check_feasibility(const typename fj_t::climber_data_t::view_t& return true; } -template -DI void save_best_l1_distance_solution(typename fj_t::climber_data_t::view_t& fj) -{ - __shared__ f_t shmem[raft::WarpSize]; - __shared__ f_t res; - // compute l1 distance on integer indices - compute_l1_distance_block( - fj.incumbent_assignment, fj.farthest_l1_sol, fj.pb.integer_indices, shmem, res); - __syncthreads(); - if (res > *fj.best_l1_distance) { - __syncthreads(); - *fj.best_l1_distance = res; - for (i_t i = threadIdx.x; i < fj.pb.n_variables; i += blockDim.x) { - fj.farthest_l1_sol[i] = fj.incumbent_assignment[i]; - } - } -} - template DI bool save_best_solution(typename fj_t::climber_data_t::view_t& fj) { cuopt_assert(blockIdx.x == 0, "Only a single block can run!"); __shared__ bool save_sol; - if (fj.settings->keep_farthest_l1_sol) { save_best_l1_distance_solution(fj); } + __shared__ bool improving; + const double improvement_threshold = max(1e-4, abs(*fj.incumbent_objective * 1e-4)); if (FIRST_THREAD) { bool better_objective = *fj.incumbent_objective < *fj.saved_solution_objective; save_sol = better_objective && fj.violated_constraints.size() == 0; + // if we're in rounding mode, save everytime we reduce the number of fractionals + if (fj.settings->mode == fj_mode_t::ROUNDING && + fj.fractional_variables.size() < *fj.saved_best_fractional_count) { + save_sol = true; + *fj.saved_best_fractional_count = fj.fractional_variables.size(); + } // save least infeasible solution if (*fj.best_excess < 0 && *fj.violation_score > *fj.best_excess) { if (fj.violated_constraints.size() == 0) @@ -615,6 +643,12 @@ DI bool save_best_solution(typename fj_t::climber_data_t::view_t& fj) *fj.best_excess = *fj.violation_score; save_sol = true; } + + if (!fj.settings->feasibility_run) { + improving = *fj.saved_solution_objective - *fj.incumbent_objective > improvement_threshold; + } else { + improving = save_sol; + } } __syncthreads(); @@ -640,10 +674,11 @@ DI bool save_best_solution(typename fj_t::climber_data_t::view_t& fj) cuopt_assert( *fj.weighted_violation_score <= *fj.max_cstr_weight * fj.pb.tolerances.absolute_tolerance, "Violated constraint and score mismatch"); - cuopt_func_call((check_feasibility(fj))); + bool check_integer = fj.settings->mode != fj_mode_t::ROUNDING; + cuopt_func_call((check_feasibility(fj, check_integer))); } // return whether it is an improving local minimum - return save_sol; + return improving; } template @@ -662,6 +697,12 @@ DI void check_exit_condition(typename fj_t::climber_data_t::view_t& vi // bool max_weight_reached = *view.max_cstr_weight > view.stop_threshold; if (local_min_cond || is_feasible) { *view.temp_break_condition = 1; } } + // if we're in rounding mode and in the "round remaining fractionals" phase, + // exit when there are no fractionals left + else if (view.settings->mode == fj_mode_t::ROUNDING && *view.handle_fractionals_only && + view.fractional_variables.size() == 0) { + *view.temp_break_condition = 1; + } } template @@ -733,14 +774,38 @@ __global__ void update_assignment_kernel(typename fj_t::climber_data_t cuopt_assert(isfinite(new_val), "assignment is not finite"); if (fj.pb.is_integer_var(var_idx)) { - cuopt_assert(fj.pb.is_integer(new_val), "The variable must be integer"); - new_val = round(new_val); + // Never "un-round" an integer variable + if (fj.pb.is_integer(fj.incumbent_assignment[var_idx])) { + cuopt_assert(fj.pb.is_integer(new_val), "The variable must be integer"); + } + + if (fj.settings->mode == fj_mode_t::ROUNDING && + !fj.pb.is_integer(fj.incumbent_assignment[var_idx]) && fj.pb.is_integer(new_val)) { + fj.fractional_variables.remove(var_idx); + + DEVICE_LOG_TRACE( + "[*][*] rounding %d from %g to %g%s, remaining %d\n", + var_idx, + fj.incumbent_assignment[var_idx], + new_val, + (new_val != round(fj.incumbent_assignment[var_idx]) ? " (non-trivial)" : " "), + fj.fractional_variables.size()); + cuopt_assert(fj.fractional_variables.size() >= 0, + "remaining fractional var count is unexpectedly negative"); + } + } + + i_t var_range = fj.pb.variable_upper_bounds[var_idx] - fj.pb.variable_lower_bounds[var_idx]; + double delta_rel_err = fabs(fj.jump_move_delta[var_idx]) / var_range; + if (delta_rel_err < fj.settings->parameters.small_move_tabu_threshold) { + *fj.small_move_tabu = *fj.iterations; } -#if defined(FJ_SINGLE_STEP) - DEVICE_LOG_TRACE( - "=---- FJ: updated %d [%g/%g] :%.4g+{%.4g}=%.4g score {%g,%g}, d_obj %.2g+%.2g=%.2g, " - "infeas %.2g, total viol %d\n", +#if FJ_SINGLE_STEP + DEVICE_LOG_DEBUG( + "=---- FJ[%d]: updated %d [%g/%g] :%.4g+{%.4g}=%.4g score {%g,%g}, d_obj %.2g+%.2g=%.2g, " + "err_range %.2g%%, infeas %.2g, total viol %d\n", + *fj.iterations, var_idx, fj.pb.variable_lower_bounds[var_idx], fj.pb.variable_upper_bounds[var_idx], @@ -752,6 +817,7 @@ __global__ void update_assignment_kernel(typename fj_t::climber_data_t *fj.incumbent_objective, fj.jump_move_delta[var_idx] * fj.pb.objective_coefficients[var_idx], *fj.incumbent_objective + fj.jump_move_delta[var_idx] * fj.pb.objective_coefficients[var_idx], + delta_rel_err, fj.jump_move_infeasibility[var_idx], fj.violated_constraints.size()); #endif @@ -825,58 +891,88 @@ DI void update_lift_moves(typename fj_t::climber_data_t::view_t fj) f_t obj_coeff = fj.pb.objective_coefficients[var_idx]; f_t delta = -std::numeric_limits::infinity(); - f_t th_lower_domain = fj.pb.variable_lower_bounds[var_idx]; - f_t th_upper_domain = fj.pb.variable_upper_bounds[var_idx]; + f_t th_lower_delta = fj.pb.variable_lower_bounds[var_idx] - fj.incumbent_assignment[var_idx]; + f_t th_upper_delta = fj.pb.variable_upper_bounds[var_idx] - fj.incumbent_assignment[var_idx]; auto [offset_begin, offset_end] = fj.pb.reverse_range_for_var(var_idx); for (i_t j = threadIdx.x + offset_begin; j < offset_end; j += blockDim.x) { - auto cstr_idx = fj.pb.reverse_constraints[j]; - auto cstr_coeff = fj.pb.reverse_coefficients[j]; - f_t c_lb = fj.pb.constraint_lower_bounds[cstr_idx]; - f_t c_ub = fj.pb.constraint_upper_bounds[cstr_idx]; + auto cstr_idx = fj.pb.reverse_constraints[j]; + auto cstr_coeff = fj.pb.reverse_coefficients[j]; + f_t c_lb = fj.pb.constraint_lower_bounds[cstr_idx]; + f_t c_ub = fj.pb.constraint_upper_bounds[cstr_idx]; + f_t cstr_tolerance = fj.get_corrected_tolerance(cstr_idx); cuopt_assert(c_lb <= c_ub, "invalid bounds"); cuopt_assert(fj.cstr_satisfied(cstr_idx, fj.incumbent_lhs[cstr_idx]), "cstr should be satisfied"); - auto [delta, sign, slack, cstr_tolerance] = - get_mtm_for_constraint( - fj, var_idx, cstr_idx, cstr_coeff); + // Process each bound separately, as both are satified and may both be finite + // otherwise range constraints aren't correctly handled + for (auto [bound, sign] : {std::make_tuple(c_lb, -1), std::make_tuple(c_ub, 1)}) { + auto [delta, slack] = + get_mtm_for_bound(fj, var_idx, cstr_idx, cstr_coeff, bound, sign); - // skip this variable if there is no slack - if (fabs(slack) <= cstr_tolerance) { - th_lower_domain = th_upper_domain = fj.incumbent_assignment[var_idx]; - continue; - } + if (cstr_coeff * sign < 0) { + if (fj.pb.is_integer_var(var_idx)) delta = ceil(delta); + } else { + if (fj.pb.is_integer_var(var_idx)) delta = floor(delta); + } - cuopt_assert(isfinite(delta), ""); - cuopt_assert(isfinite(fj.incumbent_assignment[var_idx] + delta), ""); - if (cstr_coeff * sign < 0) { - if (fj.pb.is_integer_var(var_idx)) delta = ceil(delta); - th_lower_domain = max(th_lower_domain, fj.incumbent_assignment[var_idx] + delta); - } else { - if (fj.pb.is_integer_var(var_idx)) delta = floor(delta); - th_upper_domain = min(th_upper_domain, fj.incumbent_assignment[var_idx] + delta); + // skip this variable if there is no slack + if (fabs(slack) <= cstr_tolerance) { + if (cstr_coeff * sign > 0) { + th_upper_delta = 0; + } else { + th_lower_delta = 0; + } + } else if (!fj.pb.check_variable_within_bounds(var_idx, + fj.incumbent_assignment[var_idx] + delta)) { + continue; + } else { + if (cstr_coeff * sign < 0) { + th_lower_delta = max(th_lower_delta, delta); + } else { + th_upper_delta = min(th_upper_delta, delta); + } + } } + if (th_lower_delta >= th_upper_delta) break; } // cub::BlockReduce because raft::blockReduce has a bug when using min() w/ floats // lfd = lift feasible domain - f_t lfd_lb = BlockReduce(shmem.cub).Reduce(th_lower_domain, cuda::maximum()); + f_t lfd_lb = BlockReduce(shmem.cub).Reduce(th_lower_delta, cuda::maximum()); __syncthreads(); - f_t lfd_ub = BlockReduce(shmem.cub).Reduce(th_upper_domain, cuda::minimum()); + f_t lfd_ub = BlockReduce(shmem.cub).Reduce(th_upper_delta, cuda::minimum()); - if (lfd_lb >= lfd_ub) continue; + // invalid crossing bounds + if (lfd_lb >= lfd_ub) { lfd_lb = lfd_ub = 0; } + + if (!fj.pb.check_variable_within_bounds(var_idx, fj.incumbent_assignment[var_idx] + lfd_lb)) { + lfd_lb = 0; + } + if (!fj.pb.check_variable_within_bounds(var_idx, fj.incumbent_assignment[var_idx] + lfd_ub)) { + lfd_ub = 0; + } // Now that the life move domain is computed, compute the correct lift move cuopt_assert(isfinite(fj.incumbent_assignment[var_idx]), "invalid assignment value"); delta = obj_coeff < 0 ? lfd_ub : lfd_lb; - delta = delta - fj.incumbent_assignment[var_idx]; + if (!isfinite(delta)) delta = 0; + + // check that the move is actually feasible. + cuopt_func_call(__shared__ f_t shared_delta;) if (threadIdx.x == 0) + cuopt_func_call(shared_delta = delta); + cuopt_func_call(__syncthreads()); + cuopt_func_call(auto recomputed_score = + (compute_new_score(fj, var_idx, shared_delta))); + if (threadIdx.x == 0) + cuopt_assert(recomputed_score.infeasibility >= 0, "move creates infeasibility"); // get the score auto score = fj_t::move_score_t::zero(); f_t obj_score = -1 * obj_coeff * delta; // negated to turn this into a positive score - score.base = obj_score; - if (threadIdx.x == 0 && isfinite(delta) && !fj.pb.integer_equal(delta, (f_t)0)) { + score.base = round(obj_score); + if (threadIdx.x == 0 && !fj.pb.integer_equal(delta, (f_t)0)) { fj.move_delta(FJ_MOVE_LIFT, var_idx) = delta; fj.move_score(FJ_MOVE_LIFT, var_idx) = score; fj.move_last_update(FJ_MOVE_LIFT, var_idx) = *fj.iterations; @@ -1176,6 +1272,14 @@ __global__ void select_variable_kernel(typename fj_t::climber_data_t:: auto var_idx = fj.candidate_variables.contents[setidx]; auto move_score = fj.jump_move_scores[var_idx]; + i_t var_range = fj.pb.variable_upper_bounds[var_idx] - fj.pb.variable_lower_bounds[var_idx]; + double delta_rel_err = fabs(fj.jump_move_delta[var_idx]) / var_range; + // tabu for small moves to avoid very long descents/numerical issues + if (delta_rel_err < fj.settings->parameters.small_move_tabu_threshold && + *fj.iterations - *fj.small_move_tabu < fj.settings->parameters.small_move_tabu_tenure) { + continue; + } + if (move_score > th_best_score || (move_score == th_best_score && var_idx > th_selected_var)) { th_best_score = move_score; @@ -1215,12 +1319,17 @@ __global__ void select_variable_kernel(typename fj_t::climber_data_t:: *fj.selected_var = selected_var; if (selected_var != std::numeric_limits::max()) { #if FJ_SINGLE_STEP + i_t var_range = + fj.pb.variable_upper_bounds[selected_var] - fj.pb.variable_lower_bounds[selected_var]; + double delta_rel_err = fabs(fj.jump_move_delta[selected_var]) / var_range * 100; DEVICE_LOG_INFO( - "=---- FJ: selected %d [%g/%g] :%.2g+{%.2g}=%.2g score {%g,%g}, d_obj %.2g+%.2g->%.2g, " + "=---- FJ: selected %d [%g/%g] %c :%.4g+{%.4g}=%.4g score {%g,%g}, d_obj %.2g+%.2g->%.2g, " + "delta_rel_err %.2g%%, " "infeas %.2g, total viol %d, out of %d\n", selected_var, fj.pb.variable_lower_bounds[selected_var], fj.pb.variable_upper_bounds[selected_var], + fj.pb.variable_types[selected_var] == var_t::INTEGER ? 'I' : 'C', fj.incumbent_assignment[selected_var], fj.jump_move_delta[selected_var], fj.incumbent_assignment[selected_var] + fj.jump_move_delta[selected_var], @@ -1230,6 +1339,7 @@ __global__ void select_variable_kernel(typename fj_t::climber_data_t:: fj.jump_move_delta[selected_var] * fj.pb.objective_coefficients[selected_var], *fj.incumbent_objective + fj.jump_move_delta[selected_var] * fj.pb.objective_coefficients[selected_var], + delta_rel_err, fj.jump_move_infeasibility[selected_var], fj.violated_constraints.size(), good_var_count); @@ -1384,9 +1494,12 @@ __global__ void handle_local_minimum_kernel(typename fj_t::climber_dat raft::random::PCGenerator rng(fj.settings->seed + *fj.iterations, 0, 0); __shared__ typename fj_t::move_score_t shmem[2 * raft::WarpSize]; if (*fj.break_condition) return; + // did we reach a local minimum? if (*fj.selected_var != std::numeric_limits::max()) return; + cuopt_assert(blockDim.x == TPB_localmin, "invalid TPB"); + auto best_score = fj_t::move_score_t::invalid(); i_t best_var = std::numeric_limits::max(); f_t best_delta = 0; @@ -1396,13 +1509,27 @@ __global__ void handle_local_minimum_kernel(typename fj_t::climber_dat // Local minimum, update weights. if (fj.settings->update_weights) update_weights(fj); - cuopt_assert(blockDim.x == TPB_localmin, "invalid TPB"); + // force a full refresh + *fj.full_refresh_iteration = *fj.iterations; if (blockIdx.x == 0) { bool improving = save_best_solution(fj); if (threadIdx.x == 0) { *fj.last_minimum_iteration = *fj.iterations; - if (improving) { *fj.last_improving_minimum = *fj.local_minimums_reached; } + // if we are doing a local search run (feas_run = false), + // make sure we only reset the counter when feasible solution is reached in local minimum + bool save_local_min_counter = + fj.settings->feasibility_run || fj.violated_constraints.size() == 0; + if (improving && save_local_min_counter) { + *fj.last_improving_minimum = *fj.local_minimums_reached; + // if the current objective is better than the parents (or any provided baseline) + // increase the number number of local minimax by x3 + if (fj.settings->baseline_objective_for_longer_run > + *fj.incumbent_objective + fj.settings->parameters.breakthrough_move_epsilon) { + fj.settings->n_of_minimums_for_exit *= 3; + fj.settings->baseline_objective_for_longer_run = *fj.incumbent_objective; + } + } check_exit_condition(fj); *fj.local_minimums_reached += 1; } @@ -1411,6 +1538,43 @@ __global__ void handle_local_minimum_kernel(typename fj_t::climber_dat // if we're in greedy-descent mode, stop here. if (fj.settings->mode == fj_mode_t::GREEDY_DESCENT) return; + // If we're in rounding mode and fractional variables remain, prioritize them. + if (fj.settings->mode == fj_mode_t::ROUNDING && fj.fractional_variables.size() > 0) { + if (blockIdx.x == 0) { + i_t selected = + fj.fractional_variables.contents[rng.next_u32() % fj.fractional_variables.size()]; + cuopt_assert( + fj.pb.is_integer_var(selected) && !fj.pb.is_integer(fj.incumbent_assignment[selected]), + "fractional selected variable isn't actually fractional"); + + auto [best_val, score_info] = + compute_best_mtm(fj, selected); + auto delta = best_val - fj.incumbent_assignment[selected]; + + // if no move was found, fallback to round-nearest + if (fj.pb.integer_equal(delta, 0)) { + delta = round_nearest(fj.incumbent_assignment[selected], + fj.pb.variable_lower_bounds[selected], + fj.pb.variable_upper_bounds[selected], + fj.pb.tolerances.integrality_tolerance, + rng) - + fj.incumbent_assignment[selected]; + } + + if (FIRST_THREAD) { + fj.jump_move_delta[selected] = delta; + *fj.selected_var = selected; + DEVICE_LOG_TRACE("selected_var: %d bounds [%.4g/%.4g], delta %g, old val %g\n", + *fj.selected_var, + fj.pb.variable_lower_bounds[*fj.selected_var], + fj.pb.variable_upper_bounds[*fj.selected_var], + fj.jump_move_delta[*fj.selected_var], + fj.incumbent_assignment[*fj.selected_var]); + } + } + return; + } + // Pick the best move among the variables involved in a random violated constraint. if (!fj.violated_constraints.empty()) { cg::this_grid().sync(); @@ -1470,9 +1634,6 @@ __global__ void handle_local_minimum_kernel(typename fj_t::climber_dat "assignment not within bounds"); fj.jump_move_delta[best_var] = best_delta; } - - // force a full refresh - *fj.full_refresh_iteration = *fj.iterations; } } diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cuh b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cuh index 8c839de2b..2b299ccd2 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cuh +++ b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cuh @@ -96,7 +96,7 @@ __global__ void update_changed_constraints_kernel( template __global__ void update_best_solution_kernel(typename fj_t::climber_data_t::view_t fj); -enum MTMMoveType { FJ_MTM_VIOLATED, FJ_MTM_SATISFIED }; +enum MTMMoveType { FJ_MTM_VIOLATED, FJ_MTM_SATISFIED, FJ_MTM_ALL }; template ::linear_project_onto_polytope(solution_tget_stream()); auto h_variable_lower_bounds = cuopt::host_copy(solution.problem_ptr->variable_lower_bounds, solution.handle_ptr->get_stream()); - - const f_t int_tol = context.settings.tolerances.integrality_tolerance; + auto h_last_projection = cuopt::host_copy(last_projection, solution.handle_ptr->get_stream()); + const f_t int_tol = context.settings.tolerances.integrality_tolerance; constraints_delta_t h_constraints; variables_delta_t h_variables; h_variables.n_vars = solution.problem_ptr->n_variables; @@ -182,7 +182,11 @@ bool feasibility_pump_t::linear_project_onto_polytope(solution_t constr_indices{var_id, i}; // d_j - x_j >= -val(x_j) std::vector constr_coeffs_1{1, -1}; @@ -217,15 +221,17 @@ bool feasibility_pump_t::linear_project_onto_polytope(solution_t::linear_project_onto_polytope(solution_t -bool feasibility_pump_t::random_round_with_fj(solution_t& solution, - timer_t& round_timer) -{ - const i_t n_tries = 200; - bool is_feasible = false; - rmm::device_uvector original_assign(solution.assignment, solution.handle_ptr->get_stream()); - rmm::device_uvector best_rounding(solution.assignment, solution.handle_ptr->get_stream()); - f_t best_fj_qual = std::numeric_limits::max(); - i_t i; - for (i = 0; i < n_tries; ++i) { - if (round_timer.check_time_limit()) { break; } - CUOPT_LOG_DEBUG("Trying random with FJ"); - is_feasible = solution.round_nearest(); - if (is_feasible) { - CUOPT_LOG_DEBUG("Feasible found after random round"); - return true; - } - // run fj_single descent - fj.settings.mode = fj_mode_t::GREEDY_DESCENT; - // fj.settings.n_of_minimums_for_exit = 1; - fj.settings.update_weights = false; - fj.settings.feasibility_run = true; - fj.settings.time_limit = round_timer.remaining_time(); - i_t original_sync_period = fj.settings.parameters.sync_period; - // iterations_per_graph is 10 - fj.settings.parameters.sync_period = 500; - is_feasible = fj.solve(solution); - fj.settings.parameters.sync_period = original_sync_period; - cuopt_assert(solution.test_number_all_integer(), "All integers must be rounded"); - if (is_feasible) { - CUOPT_LOG_DEBUG("Feasible found after random round FJ run"); - return true; - } - f_t current_qual = solution.get_quality(fj.cstr_weights, fj.objective_weight); - if (current_qual < best_fj_qual) { - CUOPT_LOG_DEBUG("Updating best quality of roundings %f", current_qual); - best_fj_qual = current_qual; - raft::copy(best_rounding.data(), - solution.assignment.data(), - solution.assignment.size(), - solution.handle_ptr->get_stream()); - } - raft::copy(solution.assignment.data(), - original_assign.data(), - solution.assignment.size(), - solution.handle_ptr->get_stream()); - } - n_fj_single_descents = i; - const i_t min_sols_to_run_fj = 1; - const i_t non_improving_local_minima_count = 200; - // if at least 5 solutions are explored, run longer fj - if (i >= min_sols_to_run_fj) { - // run longer fj on best sol - raft::copy(solution.assignment.data(), - best_rounding.data(), - solution.assignment.size(), - solution.handle_ptr->get_stream()); - rmm::device_uvector original_weights(fj.cstr_weights, solution.handle_ptr->get_stream()); - fj.settings.mode = fj_mode_t::EXIT_NON_IMPROVING; - fj.settings.update_weights = true; - fj.settings.feasibility_run = true; - fj.settings.time_limit = 0.5; - fj.settings.n_of_minimums_for_exit = non_improving_local_minima_count; - is_feasible = fj.solve(solution); - // restore the weights - raft::copy(fj.cstr_weights.data(), - original_weights.data(), - fj.cstr_weights.size(), - solution.handle_ptr->get_stream()); - } - if (is_feasible) { - CUOPT_LOG_DEBUG("Feasible found after %d minima FJ run", non_improving_local_minima_count); - return true; - } - raft::copy(solution.assignment.data(), - original_assign.data(), - solution.assignment.size(), - solution.handle_ptr->get_stream()); - solution.handle_ptr->sync_stream(); - return is_feasible; -} - -template -bool feasibility_pump_t::round_multiple_points(solution_t& solution) -{ - n_fj_single_descents = 0; - const f_t max_time_limit = last_lp_time * 0.1; - timer_t round_timer{std::min(max_time_limit, timer.remaining_time())}; - bool is_feasible = random_round_with_fj(solution, round_timer); - if (is_feasible) { - CUOPT_LOG_DEBUG("Feasible found after random round with fj"); - return true; - } - timer_t line_segment_timer{std::min(1., timer.remaining_time())}; - i_t n_points_to_search = n_fj_single_descents; - bool is_feasibility_run = true; - // create a copy, because assignment is changing within kernel and we want a separate point_1 - rmm::device_uvector starting_point(solution.assignment, solution.handle_ptr->get_stream()); - is_feasible = line_segment_search.search_line_segment(solution, - starting_point, - lp_optimal_solution, - n_points_to_search, - is_feasibility_run, - line_segment_timer); - if (is_feasible) { - CUOPT_LOG_DEBUG("Feasible found after line segment"); - return true; - } - // lns.config.run_lp_and_loop = false; - // timer_t lns_timer(min(last_lp_time * 0.1, timer.remaining_time())); - // is_feasible = lns.do_lns(solution, lns_timer); - // if (is_feasible) { - // CUOPT_LOG_DEBUG("Feasible found after inevitable feasibility"); - // return true; - // } - // TODO add the solution with the min distance to the population, if population is given - is_feasible = solution.round_nearest(); - if (is_feasible) { - CUOPT_LOG_DEBUG("Feasible found after nearest rounding"); - return true; - } - return is_feasible; -} - // round will use inevitable infeasibility while propagating the bounds template bool feasibility_pump_t::round(solution_t& solution) @@ -381,7 +262,7 @@ bool feasibility_pump_t::round(solution_t& solution) CUOPT_LOG_DEBUG("Rounding the point"); timer_t bounds_prop_timer(std::min(2., timer.remaining_time())); const f_t lp_run_time_after_feasible = std::min(3., timer.remaining_time() / 20.); - result = lb_constraint_prop.apply_round(solution, lp_run_time_after_feasible, bounds_prop_timer); + result = constraint_prop.apply_round(solution, lp_run_time_after_feasible, bounds_prop_timer); cuopt_func_call(solution.test_variable_bounds(true)); // copy the last rounding raft::copy(last_rounding.data(), @@ -490,7 +371,7 @@ bool feasibility_pump_t::restart_fp(solution_t& solution) fj.cstr_weights.begin(), [] __device__(f_t val) { constexpr f_t weight_divisor = 10.; - return max(f_t(10.), std::round(val / weight_divisor)); + return std::max(f_t(10.), std::round(val / weight_divisor)); }); return is_feasible; } @@ -661,15 +542,16 @@ bool feasibility_pump_t::run_single_fp_descent(solution_t& s // if the solution is almost on polytope else if (last_distances[0] < distance_to_check_for_feasible) { // run the linear projection with full precision to check if it actually is feasible - const bool return_first_feasible = true; - const f_t lp_verify_time_limit = 5.; + const f_t lp_verify_time_limit = 5.; + relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = lp_verify_time_limit; + lp_settings.tolerance = solution.problem_ptr->tolerances.absolute_tolerance; + lp_settings.return_first_feasible = true; + lp_settings.save_state = true; run_lp_with_vars_fixed(*solution.problem_ptr, solution, solution.problem_ptr->integer_indices, - context.settings.get_tolerances(), - context.lp_state, - lp_verify_time_limit, - return_first_feasible, + lp_settings, &constraint_prop.bounds_update); is_feasible = solution.get_feasible(); n_integers = solution.compute_number_of_integers(); diff --git a/cpp/src/mip/local_search/line_segment_search/line_segment_search.cu b/cpp/src/mip/local_search/line_segment_search/line_segment_search.cu index 0d1147a90..154723b3f 100644 --- a/cpp/src/mip/local_search/line_segment_search/line_segment_search.cu +++ b/cpp/src/mip/local_search/line_segment_search/line_segment_search.cu @@ -22,11 +22,13 @@ #include #include #include - +#include namespace cuopt::linear_programming::detail { template -line_segment_search_t::line_segment_search_t(fj_t& fj_) : fj(fj_) +line_segment_search_t::line_segment_search_t( + fj_t& fj_, constraint_prop_t& constraint_prop_) + : fj(fj_), constraint_prop(constraint_prop_) { } @@ -47,55 +49,155 @@ void test_point_is_within_bounds(solution_t& solution, solution.handle_ptr->get_stream()); } +class middle_first_iterator_t { + public: + middle_first_iterator_t(int n) + { + if (n > 0) queue.push({0, n - 1}); + } + + // Returns true if there is a next index + bool next(int& idx) + { + while (!queue.empty()) { + auto range = queue.front(); + queue.pop(); + int left = range.first; + int right = range.second; + if (left > right) continue; + int mid = left + (right - left) / 2; + idx = mid; + queue.push({left, mid - 1}); + queue.push({mid + 1, right}); + return true; + } + return false; + } + + private: + std::queue> queue; +}; + template -bool line_segment_search_t::search_line_segment(solution_t& solution, - const rmm::device_uvector& point_1, - const rmm::device_uvector& point_2, - i_t n_points_to_search, - bool is_feasibility_run, - cuopt::timer_t& timer) +void line_segment_search_t::save_solution_if_better( + solution_t& solution, + const rmm::device_uvector& point_1, + const rmm::device_uvector& point_2, + rmm::device_uvector& best_assignment, + rmm::device_uvector& best_feasible_assignment, + f_t& best_cost, + f_t& best_feasible_cost, + f_t curr_cost) { - CUOPT_LOG_DEBUG("Running line segment search"); + if (curr_cost < best_cost) { + bool save_solution = true; + // don't check if it is the same as parents if it is better than the parents + if (settings.recombiner_mode) { + bool is_feasible_and_better_than_parents = + solution.get_feasible() && curr_cost < settings.best_of_parents_cost; + is_feasible_and_better_than_parents = + is_feasible_and_better_than_parents || + (solution.get_feasible() && settings.parents_infeasible); + if (!is_feasible_and_better_than_parents) { + // check if the integer part is the same as one of the parents + save_solution = save_solution && !check_integer_equal_on_indices( + solution.problem_ptr->integer_indices, + solution.assignment, + point_1, + solution.problem_ptr->tolerances.integrality_tolerance, + solution.handle_ptr); + save_solution = save_solution && !check_integer_equal_on_indices( + solution.problem_ptr->integer_indices, + solution.assignment, + point_2, + solution.problem_ptr->tolerances.integrality_tolerance, + solution.handle_ptr); + } + } + if (save_solution) { + best_cost = curr_cost; + raft::copy(best_assignment.data(), + solution.assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + if (solution.get_feasible()) { + best_feasible_cost = curr_cost; + raft::copy(best_feasible_assignment.data(), + solution.assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + } + } + } +} + +template +bool line_segment_search_t::search_line_segment( + solution_t& solution, + const rmm::device_uvector& point_1, + const rmm::device_uvector& point_2, + i_t n_points_to_search, + const rmm::device_uvector& delta_vector, + bool is_feasibility_run, + cuopt::timer_t& timer) +{ + CUOPT_LOG_DEBUG("Running line segment search with a given delta vector"); cuopt_assert(point_1.size() == point_2.size(), "size mismatch"); cuopt_assert(point_1.size() == solution.assignment.size(), "size mismatch"); - cuopt_func_call(test_point_is_within_bounds(solution, point_1)); - cuopt_func_call(test_point_is_within_bounds(solution, point_2)); - rmm::device_uvector delta_vector(solution.problem_ptr->n_variables, - solution.handle_ptr->get_stream()); + cuopt_assert(delta_vector.size() == solution.assignment.size(), "size mismatch"); rmm::device_uvector best_assignment(solution.assignment, solution.handle_ptr->get_stream()); + rmm::device_uvector best_feasible_assignment(solution.assignment, + solution.handle_ptr->get_stream()); rmm::device_uvector previous_rounding(solution.assignment, solution.handle_ptr->get_stream()); - - thrust::transform(solution.handle_ptr->get_thrust_policy(), - point_2.data(), - point_2.data() + solution.problem_ptr->n_variables, - point_1.data(), - delta_vector.begin(), - [] __device__(const f_t& a, const f_t& b) { return a - b; }); - - thrust::transform( - solution.handle_ptr->get_thrust_policy(), - delta_vector.begin(), - delta_vector.end(), - delta_vector.begin(), - [n_points_to_search] __device__(const f_t& x) { return x / (n_points_to_search + 1); }); - - f_t best_cost = solution.get_quality(fj.cstr_weights, fj.objective_weight); + // we want to allow solutions that might be worse but different than parents in recombiner mode + f_t best_cost = std::numeric_limits::max(); + if (!settings.recombiner_mode) { + best_cost = solution.get_quality(fj.cstr_weights, fj.objective_weight); + } + f_t best_feasible_cost = std::numeric_limits::max(); bool initial_is_feasible = solution.get_feasible(); - // TODO start from middle and increase the resolution later - for (i_t i = 1; i <= n_points_to_search; ++i) { - CUOPT_LOG_TRACE("Line segment point %d", i); + middle_first_iterator_t it(n_points_to_search); + int i; + while (it.next(i)) { + // make it one indexed + i++; + CUOPT_LOG_DEBUG("Line segment point %d", i); thrust::tabulate(solution.handle_ptr->get_thrust_policy(), solution.assignment.begin(), solution.assignment.end(), [i, delta_ptr = delta_vector.data(), point_1_ptr = point_1.data()] __device__( const i_t index) { return point_1_ptr[index] + delta_ptr[index] * i; }); cuopt_func_call(solution.test_variable_bounds(false)); - bool is_feasible = solution.round_nearest(); - if (is_feasibility_run && is_feasible) { - CUOPT_LOG_DEBUG("Feasible found after line segment"); - return true; + bool is_feasible = false; + // if (!settings.recombiner_mode) { + if (true) { + is_feasible = solution.round_nearest(); + } else { + fj.settings.mode = fj_mode_t::ROUNDING; + fj.settings.update_weights = false; + fj.settings.feasibility_run = is_feasibility_run; + fj.settings.time_limit = std::min(0.1, timer.remaining_time()); + is_feasible = fj.solve(solution); + } + cuopt_func_call(solution.test_number_all_integer()); + if (is_feasibility_run) { + if (is_feasible) { + CUOPT_LOG_DEBUG("Feasible found during line segment round"); + return true; + } + } else { + f_t curr_cost = solution.get_quality(fj.cstr_weights, fj.objective_weight); + save_solution_if_better(solution, + point_1, + point_2, + best_assignment, + best_feasible_assignment, + best_cost, + best_feasible_cost, + curr_cost); } + if (timer.check_time_limit()) { break; } i_t number_of_integer_var_diff = compute_number_of_integer_var_diff( solution.problem_ptr->integer_indices, solution.assignment, @@ -110,10 +212,14 @@ bool line_segment_search_t::search_line_segment(solution_t& if (number_of_integer_var_diff <= min_n_integer_diffs) { continue; } cuopt_func_call(solution.test_variable_bounds(false)); // do the search here - fj.settings.mode = fj_mode_t::GREEDY_DESCENT; + fj.settings.mode = fj_mode_t::EXIT_NON_IMPROVING; + fj.settings.termination = fj_termination_flags_t::FJ_TERMINATION_ITERATION_LIMIT; + fj.settings.n_of_minimums_for_exit = 50; + fj.settings.iteration_limit = + std::max(20 * fj.settings.n_of_minimums_for_exit, solution.problem_ptr->n_constraints / 50); fj.settings.update_weights = false; fj.settings.feasibility_run = is_feasibility_run; - fj.settings.time_limit = std::min(0.5, timer.remaining_time()); + fj.settings.time_limit = std::min(1., timer.remaining_time()); is_feasible = fj.solve(solution); if (is_feasibility_run) { if (is_feasible) { @@ -122,23 +228,82 @@ bool line_segment_search_t::search_line_segment(solution_t& } } else { f_t curr_cost = solution.get_quality(fj.cstr_weights, fj.objective_weight); - if (curr_cost < best_cost && !(initial_is_feasible && !solution.get_feasible())) { - best_cost = curr_cost; - raft::copy(best_assignment.data(), - solution.assignment.data(), - solution.assignment.size(), - solution.handle_ptr->get_stream()); - } + save_solution_if_better(solution, + point_1, + point_2, + best_assignment, + best_feasible_assignment, + best_cost, + best_feasible_cost, + curr_cost); } if (timer.check_time_limit()) { break; } } - raft::copy(solution.assignment.data(), - best_assignment.data(), - solution.assignment.size(), - solution.handle_ptr->get_stream()); + // if not recombiner mode but local search mode + if (!settings.recombiner_mode) { + if (initial_is_feasible && best_feasible_cost != std::numeric_limits::max()) { + raft::copy(solution.assignment.data(), + best_feasible_assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + } else { + raft::copy(solution.assignment.data(), + best_assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + } + } else { + if (best_feasible_cost != std::numeric_limits::max()) { + CUOPT_LOG_DEBUG("Returning best feasible solution"); + // return best feasible solution that is different than parents + raft::copy(solution.assignment.data(), + best_feasible_assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + } else { + CUOPT_LOG_DEBUG("Returning best solution"); + raft::copy(solution.assignment.data(), + best_assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + } + } return solution.compute_feasibility(); } +template +bool line_segment_search_t::search_line_segment(solution_t& solution, + const rmm::device_uvector& point_1, + const rmm::device_uvector& point_2, + i_t n_points_to_search, + bool is_feasibility_run, + cuopt::timer_t& timer) +{ + CUOPT_LOG_DEBUG("Running line segment search"); + cuopt_assert(point_1.size() == point_2.size(), "size mismatch"); + cuopt_assert(point_1.size() == solution.assignment.size(), "size mismatch"); + cuopt_func_call(test_point_is_within_bounds(solution, point_1)); + cuopt_func_call(test_point_is_within_bounds(solution, point_2)); + rmm::device_uvector delta_vector(solution.problem_ptr->n_variables, + solution.handle_ptr->get_stream()); + + thrust::transform(solution.handle_ptr->get_thrust_policy(), + point_2.data(), + point_2.data() + solution.problem_ptr->n_variables, + point_1.data(), + delta_vector.begin(), + [] __device__(const f_t a, const f_t b) { return a - b; }); + + thrust::transform( + solution.handle_ptr->get_thrust_policy(), + delta_vector.begin(), + delta_vector.end(), + delta_vector.begin(), + [n_points_to_search] __device__(const f_t x) { return x / (n_points_to_search + 1); }); + return search_line_segment( + solution, point_1, point_2, n_points_to_search, delta_vector, is_feasibility_run, timer); +} + #if MIP_INSTANTIATE_FLOAT template class line_segment_search_t; #endif diff --git a/cpp/src/mip/local_search/line_segment_search/line_segment_search.cuh b/cpp/src/mip/local_search/line_segment_search/line_segment_search.cuh index 70a0b9193..a7ee6ef15 100644 --- a/cpp/src/mip/local_search/line_segment_search/line_segment_search.cuh +++ b/cpp/src/mip/local_search/line_segment_search/line_segment_search.cuh @@ -18,23 +18,49 @@ #pragma once #include +#include #include namespace cuopt::linear_programming::detail { +struct line_segment_settings_t { + bool recombiner_mode = false; + double best_of_parents_cost = std::numeric_limits::max(); + bool parents_infeasible = false; +}; + template class line_segment_search_t { public: line_segment_search_t() = delete; - line_segment_search_t(fj_t& fj); + line_segment_search_t(fj_t& fj, constraint_prop_t& constraint_prop); + bool search_line_segment(solution_t& solution, + const rmm::device_uvector& point_1, + const rmm::device_uvector& point_2, + i_t n_points_to_search, + bool is_feasibility_run, + cuopt::timer_t& timer); + bool search_line_segment(solution_t& solution, const rmm::device_uvector& point_1, const rmm::device_uvector& point_2, i_t n_points_to_search, + const rmm::device_uvector& delta_vector, bool is_feasibility_run, cuopt::timer_t& timer); + void save_solution_if_better(solution_t& solution, + const rmm::device_uvector& point_1, + const rmm::device_uvector& point_2, + rmm::device_uvector& best_assignment, + rmm::device_uvector& best_feasible_assignment, + f_t& best_cost, + f_t& best_feasible_cost, + f_t curr_cost); + fj_t& fj; + constraint_prop_t& constraint_prop; + line_segment_settings_t settings; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/local_search/local_search.cu b/cpp/src/mip/local_search/local_search.cu index 7c02fefdc..655b07fa9 100644 --- a/cpp/src/mip/local_search/local_search.cu +++ b/cpp/src/mip/local_search/local_search.cu @@ -41,7 +41,7 @@ local_search_t::local_search_t(mip_solver_context_t& context // fj_tree(fj), constraint_prop(context), lb_constraint_prop(context), - line_segment_search(fj), + line_segment_search(fj, constraint_prop), fp(context, fj, // fj_tree, @@ -85,38 +85,36 @@ void local_search_t::generate_fast_solution(solution_t& solu template bool local_search_t::run_local_search(solution_t& solution, const weight_t& weights, - timer_t timer) + timer_t timer, + f_t baseline_objective, + bool at_least_one_parent_feasible) { raft::common::nvtx::range fun_scope("local search"); fj_settings_t fj_settings; if (timer.check_time_limit()) return false; // adjust these time limits - fj_settings.time_limit = timer.remaining_time(); + if (!solution.get_feasible()) { + if (at_least_one_parent_feasible) { + fj_settings.time_limit = 1.; + timer = timer_t(1.); + } else { + fj_settings.time_limit = 0.5; + timer = timer_t(0.5); + } + } else { + fj_settings.time_limit = timer.remaining_time(); + } fj_settings.update_weights = false; fj_settings.feasibility_run = false; fj.set_fj_settings(fj_settings); fj.copy_weights(weights, solution.handle_ptr); - cuopt_func_call(bool is_feasible_before_search = solution.get_feasible()); bool is_feas; i_t rd = std::uniform_int_distribution(0, 1)(rng); if (rd == 0 && lp_optimal_exists) { is_feas = run_fj_line_segment(solution, timer); - cuopt_assert(!is_feasible_before_search || is_feas, - "Line segment search should not change feasibility"); } else { - solution_t solution_before_annealing(solution); - bool is_feasible_before_annealing = solution.get_feasible(); - is_feas = run_fj_annealing(solution, timer); - // we should remove this logic once, FJ is fixed - if (is_feasible_before_annealing && !solution.get_feasible()) { - solution.copy_from(solution_before_annealing); - solution.handle_ptr->sync_stream(); - } - cuopt_func_call(bool is_feasible_after_annealing = solution.get_feasible()); - cuopt_assert(!is_feasible_before_search || is_feasible_after_annealing, - "Annealing should not change feasibility"); + is_feas = run_fj_annealing(solution, timer, baseline_objective); } - // TODO add a test for the quality improvement too return is_feas; } @@ -141,18 +139,23 @@ bool local_search_t::run_fj_until_timer(solution_t& solution // SIMULATED ANNEALING not fully implemented yet, placeholder template -bool local_search_t::run_fj_annealing(solution_t& solution, timer_t timer) +bool local_search_t::run_fj_annealing(solution_t& solution, + timer_t timer, + f_t baseline_objective) { auto prev_settings = fj.settings; // run in FEASIBLE_FIRST to priorize feasibility-improving moves - fj.settings.n_of_minimums_for_exit = 50; + fj.settings.n_of_minimums_for_exit = 250; fj.settings.mode = fj_mode_t::EXIT_NON_IMPROVING; fj.settings.candidate_selection = fj_candidate_selection_t::FEASIBLE_FIRST; - fj.settings.termination = fj_termination_flags_t::FJ_TERMINATION_TIME_LIMIT; - fj.settings.time_limit = std::min(10., timer.remaining_time()); + fj.settings.termination = fj_termination_flags_t::FJ_TERMINATION_ITERATION_LIMIT; + fj.settings.iteration_limit = + max(20 * fj.settings.n_of_minimums_for_exit, solution.problem_ptr->n_constraints / 50); + fj.settings.time_limit = std::min(10., timer.remaining_time()); fj.settings.parameters.allow_infeasibility_iterations = 100; fj.settings.update_weights = 1; + fj.settings.baseline_objective_for_longer_run = baseline_objective; fj.solve(solution); bool is_feasible = solution.compute_feasibility(); @@ -194,12 +197,11 @@ bool local_search_t::check_fj_on_lp_optimal(solution_t& solu constraint_prop.apply_round(solution, lp_run_time_after_feasible, bounds_prop_timer); if (!is_feasible) { const f_t lp_run_time = 2.; - run_lp_with_vars_fixed(*solution.problem_ptr, - solution, - solution.problem_ptr->integer_indices, - solution.problem_ptr->tolerances, - context.lp_state, - lp_run_time); + relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = lp_run_time; + lp_settings.tolerance = solution.problem_ptr->tolerances.absolute_tolerance; + run_lp_with_vars_fixed( + *solution.problem_ptr, solution, solution.problem_ptr->integer_indices, lp_settings); } else { return is_feasible; } diff --git a/cpp/src/mip/local_search/local_search.cuh b/cpp/src/mip/local_search/local_search.cuh index 9eeeb37ee..4334f1063 100644 --- a/cpp/src/mip/local_search/local_search.cuh +++ b/cpp/src/mip/local_search/local_search.cuh @@ -41,8 +41,12 @@ class local_search_t { timer_t timer); bool run_local_search(solution_t& solution, const weight_t& weights, - timer_t timer); - bool run_fj_annealing(solution_t& solution, timer_t timer); + timer_t timer, + f_t baseline_objective = std::numeric_limits::lowest(), + bool at_least_one_parent_feasible = true); + bool run_fj_annealing(solution_t& solution, + timer_t timer, + f_t baseline_objective = std::numeric_limits::lowest()); bool run_fj_line_segment(solution_t& solution, timer_t timer); bool run_fj_on_zero(solution_t& solution, timer_t timer); bool check_fj_on_lp_optimal(solution_t& solution, bool perturb, timer_t timer); diff --git a/cpp/src/mip/local_search/local_search_config.hpp b/cpp/src/mip/local_search/local_search_config.hpp new file mode 100644 index 000000000..56e17c30a --- /dev/null +++ b/cpp/src/mip/local_search/local_search_config.hpp @@ -0,0 +1,29 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +namespace cuopt::linear_programming::detail { + +struct ls_config_t { + static constexpr bool use_line_segment = true; + static constexpr bool use_fj = true; + static constexpr bool use_fp_ls = false; + static constexpr bool use_cutting_plane_from_best_solution = false; +}; + +} // namespace cuopt::linear_programming::detail \ No newline at end of file diff --git a/cpp/src/mip/local_search/rounding/constraint_prop.cu b/cpp/src/mip/local_search/rounding/constraint_prop.cu index e681474ad..33a631d29 100644 --- a/cpp/src/mip/local_search/rounding/constraint_prop.cu +++ b/cpp/src/mip/local_search/rounding/constraint_prop.cu @@ -39,6 +39,7 @@ constraint_prop_t::constraint_prop_t(mip_solver_context_t& c temp_problem(*context.problem_ptr), temp_sol(*context.problem_ptr), bounds_update(context), + multi_probe(context), bounds_repair(*context.problem_ptr, bounds_update), conditional_bounds_update(*context.problem_ptr), set_vars(context.problem_ptr->n_variables, context.problem_ptr->handle_ptr->get_stream()), @@ -176,12 +177,16 @@ void constraint_prop_t::sort_by_implied_slack_consumption(solution_t implied_slack_consumption_per_var(vars.size(), sol.handle_ptr->get_stream()); const i_t block_dim = 128; + auto min_activity = selected_update ? make_span(multi_probe.upd_1.min_activity) + : make_span(multi_probe.upd_0.min_activity); + auto max_activity = selected_update ? make_span(multi_probe.upd_1.max_activity) + : make_span(multi_probe.upd_0.max_activity); compute_implied_slack_consumption_per_var <<get_stream()>>>( sol.problem_ptr->view(), vars, - make_span(bounds_update.upd.min_activity), - make_span(bounds_update.upd.max_activity), + min_activity, + max_activity, make_span(implied_slack_consumption_per_var), problem_ii, context.settings.get_tolerances()); @@ -516,14 +521,24 @@ struct greater_than_threshold_t { template void constraint_prop_t::copy_bounds(rmm::device_uvector& output_lb, rmm::device_uvector& output_ub, - rmm::device_uvector& output_assignment, const rmm::device_uvector& input_lb, const rmm::device_uvector& input_ub, - const rmm::device_uvector& input_assignment, const raft::handle_t* handle_ptr) { raft::copy(output_lb.data(), input_lb.data(), input_lb.size(), handle_ptr->get_stream()); raft::copy(output_ub.data(), input_ub.data(), input_ub.size(), handle_ptr->get_stream()); +} + +template +void constraint_prop_t::copy_bounds(rmm::device_uvector& output_lb, + rmm::device_uvector& output_ub, + rmm::device_uvector& output_assignment, + const rmm::device_uvector& input_lb, + const rmm::device_uvector& input_ub, + const rmm::device_uvector& input_assignment, + const raft::handle_t* handle_ptr) +{ + copy_bounds(output_lb, output_ub, input_lb, input_ub, handle_ptr); raft::copy(output_assignment.data(), input_assignment.data(), input_assignment.size(), @@ -568,25 +583,83 @@ void constraint_prop_t::restore_original_bounds(solution_t& } template -std::vector> constraint_prop_t::generate_bulk_rounding_vector( +thrust::pair constraint_prop_t::generate_double_probing_pair( + const solution_t& sol, + const solution_t& orig_sol, + i_t unset_var_idx, + const std::optional>> probing_config, + bool bulk_rounding) +{ + f_t first_probe, second_probe; + if (probing_config.has_value()) { + // for now get the first one + auto [from_first, from_second] = probing_config.value().get().probing_values[unset_var_idx]; + std::mt19937 rng(cuopt::seed_generator::get_seed()); + std::uniform_real_distribution dist(0.0f, 1.0f); + f_t random_value = dist(rng); + f_t average_value = (from_first + from_second) / 2; + if (random_value > 0.5) { + average_value = ceil(average_value); + } else { + average_value = floor(average_value); + } + + // if the first set of values are more represented use the second value + if (probing_config.value().get().n_of_fixed_from_first > + probing_config.value().get().n_of_fixed_from_second) { + // if it is the same as less represented + if (average_value == from_second) { + first_probe = from_second; + second_probe = from_first; + } else { + first_probe = average_value; + second_probe = from_second; + } + } else { + // if it is the same as less represented + if (average_value == from_first) { + first_probe = from_first; + second_probe = from_second; + } else { + first_probe = average_value; + second_probe = from_first; + } + } + } else { + std::tie(first_probe, std::ignore, second_probe) = probing_values(sol, orig_sol, unset_var_idx); + // do another draw in bulk rounding, so the second vector is randomly drawn + if (bulk_rounding) { + std::tie(second_probe, std::ignore, std::ignore) = + probing_values(sol, orig_sol, unset_var_idx); + } + } + return thrust::make_pair(first_probe, second_probe); +} + +template +std::tuple, std::vector, std::vector> +constraint_prop_t::generate_bulk_rounding_vector( const solution_t& sol, const solution_t& orig_sol, const std::vector& host_vars_to_set, - const std::optional>> probing_candidates) + const std::optional>> probing_config) { const f_t int_tol = orig_sol.problem_ptr->tolerances.integrality_tolerance; std::string log_str{"Setting var:\t"}; - std::vector> var_val_pairs; - var_val_pairs.reserve(host_vars_to_set.size()); + std::tuple, std::vector, std::vector> var_probe_vals; + std::get<0>(var_probe_vals).resize(host_vars_to_set.size()); + std::get<1>(var_probe_vals).resize(host_vars_to_set.size()); + std::get<2>(var_probe_vals).resize(host_vars_to_set.size()); for (i_t i = 0; i < (i_t)host_vars_to_set.size(); ++i) { auto unset_var_idx = host_vars_to_set[i]; f_t first_probe, second_probe; - if (probing_candidates.has_value()) { - // for now get the first one - thrust::tie(first_probe, second_probe) = probing_candidates.value()[unset_var_idx]; + // if it is a bulk rounding do + if (host_vars_to_set.size() > 1) { + cuda::std::tie(first_probe, second_probe) = + generate_double_probing_pair(sol, orig_sol, unset_var_idx, probing_config, true); } else { - std::tie(first_probe, std::ignore, second_probe) = - probing_values(sol, orig_sol, unset_var_idx); + cuda::std::tie(first_probe, second_probe) = + generate_double_probing_pair(sol, orig_sol, unset_var_idx, probing_config, false); } cuopt_assert(orig_sol.problem_ptr->is_integer(first_probe), "Probing value must be an integer"); cuopt_assert(orig_sol.problem_ptr->is_integer(second_probe), @@ -596,25 +669,32 @@ std::vector> constraint_prop_t::generate_bulk_r if (use_probing_cache && bounds_update.probing_cache.contains(*sol.problem_ptr, unset_var_idx)) { // check if there are any conflicting bounds - val_to_round = - bounds_update.probing_cache.get_least_conflicting_rounding(*sol.problem_ptr, - bounds_update.host_lb, - bounds_update.host_ub, - unset_var_idx, - first_probe, - second_probe, - int_tol); + val_to_round = bounds_update.probing_cache.get_least_conflicting_rounding(*sol.problem_ptr, + multi_probe.host_lb, + multi_probe.host_ub, + unset_var_idx, + first_probe, + second_probe, + int_tol); + if (val_to_round == second_probe) { second_probe = first_probe; } } cuopt_assert(orig_sol.problem_ptr->variable_lower_bounds.element( unset_var_idx, sol.handle_ptr->get_stream()) <= val_to_round + int_tol && val_to_round - int_tol <= orig_sol.problem_ptr->variable_upper_bounds.element( unset_var_idx, sol.handle_ptr->get_stream()), "Variable out of original bounds!"); - var_val_pairs.emplace_back(unset_var_idx, val_to_round); + cuopt_assert(orig_sol.problem_ptr->variable_lower_bounds.element( + unset_var_idx, sol.handle_ptr->get_stream()) <= second_probe + int_tol && + second_probe - int_tol <= orig_sol.problem_ptr->variable_upper_bounds.element( + unset_var_idx, sol.handle_ptr->get_stream()), + "Variable out of original bounds!"); + std::get<0>(var_probe_vals)[i] = unset_var_idx; + std::get<1>(var_probe_vals)[i] = val_to_round; + std::get<2>(var_probe_vals)[i] = second_probe; log_str.append(std::to_string(unset_var_idx) + ", "); } CUOPT_LOG_TRACE("%s", log_str.c_str()); - return var_val_pairs; + return var_probe_vals; } template @@ -629,15 +709,15 @@ void constraint_prop_t::update_host_assignment(const solution_t void constraint_prop_t::set_host_bounds(const solution_t& sol) { - cuopt_assert(sol.problem_ptr->variable_lower_bounds.size() == bounds_update.host_lb.size(), + cuopt_assert(sol.problem_ptr->variable_lower_bounds.size() == multi_probe.host_lb.size(), "size of variable lower bound mismatch"); - raft::copy(bounds_update.host_lb.data(), + raft::copy(multi_probe.host_lb.data(), sol.problem_ptr->variable_lower_bounds.data(), sol.problem_ptr->variable_lower_bounds.size(), sol.handle_ptr->get_stream()); - cuopt_assert(sol.problem_ptr->variable_upper_bounds.size() == bounds_update.host_ub.size(), + cuopt_assert(sol.problem_ptr->variable_upper_bounds.size() == multi_probe.host_ub.size(), "size of variable upper bound mismatch"); - raft::copy(bounds_update.host_ub.data(), + raft::copy(multi_probe.host_ub.data(), sol.problem_ptr->variable_upper_bounds.data(), sol.problem_ptr->variable_upper_bounds.size(), sol.handle_ptr->get_stream()); @@ -668,7 +748,10 @@ bool constraint_prop_t::run_repair_procedure(problem_t& prob timer_t& timer, const raft::handle_t* handle_ptr) { - bounds_update.set_updated_bounds(problem); + // select the first probing value + i_t select = 0; + multi_probe.set_updated_bounds(problem, select, handle_ptr); + bounds_update.copy_input_bounds(problem); repair_stats.repair_attempts++; f_t repair_start_time = timer.remaining_time(); i_t n_of_repairs_needed_for_feasible = 0; @@ -693,8 +776,9 @@ bool constraint_prop_t::run_repair_procedure(problem_t& prob // it. note that the number of fixed vars will still be the same, as repair only shifts the vars restore_original_bounds_on_unfixed(problem, original_problem, handle_ptr); bounds_update.settings.iteration_limit = 100; + bounds_update.settings.time_limit = timer.remaining_time(); auto term_crit = bounds_update.solve(problem); - bounds_update.settings.iteration_limit = 20; + bounds_update.settings.iteration_limit = 50; if (timer.check_time_limit()) { CUOPT_LOG_DEBUG("Time limit is reached in repair loop!"); f_t repair_end_time = timer.remaining_time(); @@ -713,6 +797,8 @@ bool constraint_prop_t::run_repair_procedure(problem_t& prob CUOPT_LOG_DEBUG("Repair success: n_of_repair_calls needed: %d", n_of_repairs_needed_for_feasible); f_t repair_end_time = timer.remaining_time(); repair_stats.total_time_spent_on_repair += repair_start_time - repair_end_time; + // test that bounds are really repaired and no ii cstr is present + cuopt_assert(!is_problem_ii(problem), "Problem must not be ii after repair success"); return true; } @@ -731,7 +817,7 @@ bool constraint_prop_t::find_integer( solution_t& orig_sol, f_t lp_run_time_after_feasible, timer_t& timer, - std::optional>> probing_candidates) + std::optional>> probing_config) { using crit_t = termination_criterion_t; auto& unset_integer_vars = unset_vars; @@ -743,8 +829,11 @@ bool constraint_prop_t::find_integer( curr_host_assignment.resize(sol.problem_ptr->n_variables); // round vals that are close enough bounds_update.settings.time_limit = max_timer.remaining_time(); - bounds_update.settings.iteration_limit = 20; + bounds_update.settings.iteration_limit = 50; bounds_update.resize(*sol.problem_ptr); + multi_probe.settings.iteration_limit = 50; + multi_probe.settings.time_limit = max_timer.remaining_time(); + multi_probe.resize(*sol.problem_ptr); if (max_timer.check_time_limit()) { CUOPT_LOG_DEBUG("Time limit is reached before bounds prop rounding!"); sol.round_nearest(); @@ -786,10 +875,11 @@ bool constraint_prop_t::find_integer( sort_by_interval_and_frac(sol, make_span(unset_integer_vars), rng); } set_host_bounds(sol); - size_t set_count = 0; - bool timeout_happened = false; - bool repair_tried = false; + size_t set_count = 0; + bool timeout_happened = false; + i_t n_failed_repair_iterations = 0; while (set_count < unset_integer_vars.size()) { + CUOPT_LOG_TRACE("n_set_vars %d vars to set %lu", set_count, unset_integer_vars.size()); update_host_assignment(sol); if (max_timer.check_time_limit()) { CUOPT_LOG_DEBUG("Second time limit is reached returning nearest rounding!"); @@ -800,11 +890,13 @@ bool constraint_prop_t::find_integer( if (!rounding_ii && timer.check_time_limit()) { CUOPT_LOG_DEBUG("First time limit is reached! Continuing without backtracking and repair!"); rounding_ii = true; - // this is for not trying the repair procedure again - repair_tried = true; + // this is to not try the repair procedure again + n_failed_repair_iterations = std::numeric_limits::max(); } const i_t n_curr_unset = unset_integer_vars.size() - set_count; - if (!recovery_mode || rounding_ii) { + if (single_rounding_only) { + bounds_prop_interval = 1; + } else if (!recovery_mode || rounding_ii) { if (n_curr_unset > 36) { bounds_prop_interval = sqrt(n_curr_unset); } else { @@ -822,10 +914,12 @@ bool constraint_prop_t::find_integer( unset_integer_vars.data() + set_count, n_vars_to_set, sol.handle_ptr->get_stream()); - auto var_val_pairs = - generate_bulk_rounding_vector(sol, orig_sol, host_vars_to_set, probing_candidates); - probe(sol, orig_sol.problem_ptr, var_val_pairs, &set_count, unset_integer_vars); - if (!repair_tried && rounding_ii && !timeout_happened) { + auto var_probe_vals = + generate_bulk_rounding_vector(sol, orig_sol, host_vars_to_set, probing_config); + probe( + sol, orig_sol.problem_ptr, var_probe_vals, &set_count, unset_integer_vars, probing_config); + if (!(n_failed_repair_iterations >= max_n_failed_repair_iterations) && rounding_ii && + !timeout_happened) { timer_t repair_timer{std::min(timer.remaining_time() / 5, timer.elapsed_time() / 3)}; save_bounds(sol); // update bounds and run repair procedure @@ -833,7 +927,7 @@ bool constraint_prop_t::find_integer( run_repair_procedure(*sol.problem_ptr, *orig_sol.problem_ptr, repair_timer, sol.handle_ptr); if (!bounds_repaired) { restore_bounds(sol); - repair_tried = true; + n_failed_repair_iterations++; } else { CUOPT_LOG_DEBUG( "Bounds are repaired! Deactivating recovery mode in bounds prop. n_curr_unset %d " @@ -843,9 +937,6 @@ bool constraint_prop_t::find_integer( recovery_mode = false; rounding_ii = false; n_iter_in_recovery = 0; - // test that bounds are really repaired and no ii cstr is present - cuopt_assert(!is_problem_ii(*sol.problem_ptr), - "Problem must not be ii after repair success"); // during repair procedure some variables might be collapsed auto iter = thrust::stable_partition( sol.handle_ptr->get_thrust_policy(), @@ -862,8 +953,9 @@ bool constraint_prop_t::find_integer( set_count += n_fixed_vars; } } - if (recovery_mode && bounds_update.infeas_constraints_count > 0) { - // if bounds are not repaired, restore previous bounds + // we can keep normal bounds_update here because this is only activated after the repair + if (recovery_mode && multi_probe.infeas_constraints_count_0 > 0 && + multi_probe.infeas_constraints_count_1 > 0) { CUOPT_LOG_DEBUG("Problem is ii in constraint prop. n_curr_unset %d bounds_prop_interval %d", n_curr_unset, bounds_prop_interval); @@ -882,24 +974,31 @@ bool constraint_prop_t::find_integer( // we update from the problem bounds and not the final bounds of bounds update // because we might be in a recovery mode where we want to continue with the bounds before bulk // which is the unchanged problem bounds - bounds_update.update_host_bounds(sol.handle_ptr, - make_span(sol.problem_ptr->variable_lower_bounds), - make_span(sol.problem_ptr->variable_upper_bounds)); + multi_probe.update_host_bounds(sol.handle_ptr, + make_span(sol.problem_ptr->variable_lower_bounds), + make_span(sol.problem_ptr->variable_upper_bounds)); } - CUOPT_LOG_DEBUG("Bounds propagation rounding end: ii constraint count %d", - bounds_update.infeas_constraints_count); + CUOPT_LOG_DEBUG( + "Bounds propagation rounding end: ii constraint count first buffer %d, second buffer %d", + multi_probe.infeas_constraints_count_0, + multi_probe.infeas_constraints_count_1); cuopt_assert(sol.test_number_all_integer(), "All integers must be rounded"); expand_device_copy(orig_sol.assignment, sol.assignment, sol.handle_ptr->get_stream()); cuopt_func_call(orig_sol.test_variable_bounds()); // if the constraint is not ii, run LP - if (bounds_update.infeas_constraints_count == 0 && !timeout_happened) { + if ((multi_probe.infeas_constraints_count_0 == 0 || + multi_probe.infeas_constraints_count_1 == 0) && + !timeout_happened) { + relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = lp_run_time_after_feasible; + lp_settings.tolerance = orig_sol.problem_ptr->tolerances.absolute_tolerance; + lp_settings.save_state = false; + lp_settings.return_first_feasible = true; run_lp_with_vars_fixed(*orig_sol.problem_ptr, orig_sol, orig_sol.problem_ptr->integer_indices, - context.settings.get_tolerances(), - context.lp_state, - lp_run_time_after_feasible, - true); + lp_settings, + static_cast*>(nullptr)); } bool res_feasible = orig_sol.compute_feasibility(); orig_sol.handle_ptr->sync_stream(); @@ -911,7 +1010,7 @@ bool constraint_prop_t::apply_round( solution_t& sol, f_t lp_run_time_after_feasible, timer_t& timer, - std::optional>> probing_candidates) + std::optional>> probing_config) { raft::common::nvtx::range fun_scope("constraint prop round"); @@ -928,8 +1027,7 @@ bool constraint_prop_t::apply_round( temp_sol.problem_ptr = &p; f_t bounds_prop_start_time = max_timer.remaining_time(); cuopt_func_call(temp_sol.test_variable_bounds(false)); - bool sol_found = - find_integer(temp_sol, sol, lp_run_time_after_feasible, timer, probing_candidates); + bool sol_found = find_integer(temp_sol, sol, lp_run_time_after_feasible, timer, probing_config); f_t bounds_prop_end_time = max_timer.remaining_time(); repair_stats.total_time_spent_on_bounds_prop += bounds_prop_start_time - bounds_prop_end_time; @@ -955,8 +1053,8 @@ template std::tuple constraint_prop_t::probing_values( const solution_t& sol, const solution_t& orig_sol, i_t idx) { - auto v_lb = bounds_update.host_lb[idx]; - auto v_ub = bounds_update.host_ub[idx]; + auto v_lb = multi_probe.host_lb[idx]; + auto v_ub = multi_probe.host_ub[idx]; auto var_val = curr_host_assignment[idx]; const f_t int_tol = sol.problem_ptr->tolerances.integrality_tolerance; @@ -1020,34 +1118,15 @@ std::tuple constraint_prop_t::probing_values( } template -bool constraint_prop_t::probe( +bool constraint_prop_t::handle_fixed_vars( solution_t& sol, problem_t* original_problem, - const std::vector>& var_probe_val_pairs, + const std::tuple, std::vector, std::vector>& var_probe_vals, size_t* set_count_ptr, rmm::device_uvector& unset_vars) { - const f_t int_tol = sol.problem_ptr->tolerances.integrality_tolerance; - auto set_count = *set_count_ptr; - const bool use_host_bounds = true; - bounds_update.solve(*sol.problem_ptr, var_probe_val_pairs, use_host_bounds); - // if we are ii at this point, backtrack the number of variables we have set in this given - // interval then start setting one by one - // if we determined that the rounding is ii then don't do any recovery and finish ronuding - // quickly - bool bounds_update_ii = bounds_update.infeas_constraints_count > 0; - if (!recovery_mode && !rounding_ii && bounds_update_ii && bounds_prop_interval != 1) { - CUOPT_LOG_DEBUG("Activating recovery mode in bounds prop: bounds_prop_interval %d n_ii cstr %d", - bounds_prop_interval, - bounds_update.infeas_constraints_count); - // do backtracking - recovery_mode = true; - bounds_update.infeas_constraints_count = 0; - n_iter_in_recovery = 0; - return false; - } - bounds_update.set_updated_bounds(*sol.problem_ptr); - + auto set_count = *set_count_ptr; + const f_t int_tol = sol.problem_ptr->tolerances.integrality_tolerance; // which other variables were affected? auto iter = thrust::stable_partition( sol.handle_ptr->get_thrust_policy(), @@ -1060,44 +1139,71 @@ bool constraint_prop_t::probe( make_span(original_problem->variable_upper_bounds), make_span(sol.assignment)}); i_t n_fixed_vars = (iter - (unset_vars.begin() + set_count)); - cuopt_assert(n_fixed_vars >= var_probe_val_pairs.size(), "Error in number of vars fixed!"); + cuopt_assert(n_fixed_vars >= std::get<0>(var_probe_vals).size(), + "Error in number of vars fixed!"); set_count += n_fixed_vars; CUOPT_LOG_TRACE("Set var count increased from %d to %d", *set_count_ptr, set_count); *set_count_ptr = set_count; - return bounds_update.infeas_constraints_count == 0; + return multi_probe.infeas_constraints_count_0 == 0 || multi_probe.infeas_constraints_count_1 == 0; } template -void constraint_prop_t::relax_crossing_bound_vars(solution_t& sol, - raft::device_span lower_bounds, - raft::device_span upper_bounds) +bool constraint_prop_t::probe( + solution_t& sol, + problem_t* original_problem, + const std::tuple, std::vector, std::vector>& var_probe_vals, + size_t* set_count_ptr, + rmm::device_uvector& unset_vars, + std::optional>> probing_config) { - const f_t int_tol = sol.problem_ptr->tolerances.integrality_tolerance; - - thrust::for_each(sol.handle_ptr->get_thrust_policy(), - thrust::make_counting_iterator(0), - thrust::make_counting_iterator(sol.problem_ptr->n_variables), - [ass_ptr = sol.assignment.data(), - int_tol, - lower_bounds, - upper_bounds, - var_type_ptr = sol.problem_ptr->variable_types.data()] __device__(i_t idx) { - if (lower_bounds[idx] - int_tol > upper_bounds[idx]) { - printf("relaxing var with %d with bounds %f and %f var type %d\n", - idx, - lower_bounds[idx], - upper_bounds[idx], - (i_t)var_type_ptr[idx]); - // set one of finite bounds - if (!isfinite(lower_bounds[idx])) { - ass_ptr[idx] = ceil(upper_bounds[idx] - 0.5); - } else if (!isfinite(upper_bounds[idx])) { - ass_ptr[idx] = floor(lower_bounds[idx] + 0.5); - } else { - ass_ptr[idx] = round((lower_bounds[idx] + upper_bounds[idx]) / 2); - } - } - }); + const bool use_host_bounds = true; + multi_probe.settings.time_limit = max_timer.remaining_time(); + multi_probe.solve(*sol.problem_ptr, var_probe_vals, use_host_bounds); + // if we are ii at this point, backtrack the number of variables we have set in this given + // interval then start setting one by one + // if we determined that the rounding is ii then don't do any recovery and finish ronuding + // quickly + bool first_bounds_update_ii = multi_probe.infeas_constraints_count_0 > 0; + bool second_bounds_update_ii = multi_probe.infeas_constraints_count_1 > 0; + // if we are on single rounding mode, direcly mark as ii so we can run repair + if (!rounding_ii && single_rounding_only && (first_bounds_update_ii && second_bounds_update_ii)) { + rounding_ii = true; + return false; + } + if (!recovery_mode && !rounding_ii && (first_bounds_update_ii && second_bounds_update_ii) && + bounds_prop_interval != 1) { + CUOPT_LOG_DEBUG( + "Activating recovery mode in bounds prop: bounds_prop_interval %d n_ii cstr %d %d", + bounds_prop_interval, + multi_probe.infeas_constraints_count_0, + multi_probe.infeas_constraints_count_1); + // do backtracking + recovery_mode = true; + multi_probe.infeas_constraints_count_0 = 0; + multi_probe.infeas_constraints_count_1 = 0; + n_iter_in_recovery = 0; + return false; + } + selected_update = 0; + if (first_bounds_update_ii) { selected_update = 1; } + // if we are doing single rounding + if (probing_config.has_value() && probing_config.value().get().use_balanced_probing) { + cuopt_assert(std::get<0>(var_probe_vals).size() == 1, + "Balanced probing must be used with single rounding"); + i_t var_idx = std::get<0>(var_probe_vals)[0]; + f_t value_chosen = + selected_update ? std::get<1>(var_probe_vals)[0] : std::get<2>(var_probe_vals)[0]; + if (value_chosen == probing_config.value().get().probing_values[var_idx].first) { + probing_config.value().get().n_of_fixed_from_first++; + } else { + probing_config.value().get().n_of_fixed_from_second++; + } + CUOPT_LOG_TRACE("Balanced probing: n_of_fixed_from_first %d n_of_fixed_from_second %d", + probing_config.value().get().n_of_fixed_from_first, + probing_config.value().get().n_of_fixed_from_second); + } + multi_probe.set_updated_bounds(*sol.problem_ptr, selected_update, sol.handle_ptr); + return handle_fixed_vars(sol, original_problem, var_probe_vals, set_count_ptr, unset_vars); } #if MIP_INSTANTIATE_FLOAT diff --git a/cpp/src/mip/local_search/rounding/constraint_prop.cuh b/cpp/src/mip/local_search/rounding/constraint_prop.cuh index 1adea878d..cb9cbbd00 100644 --- a/cpp/src/mip/local_search/rounding/constraint_prop.cuh +++ b/cpp/src/mip/local_search/rounding/constraint_prop.cuh @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -37,14 +38,23 @@ struct repair_stats_t { double total_time_spent_on_bounds_prop = 0.; }; +template +struct probing_config_t { + probing_config_t(i_t n_vars, const raft::handle_t* handle_ptr) : probing_values(n_vars) {} + bool use_balanced_probing = false; + i_t n_of_fixed_from_first = 0; + i_t n_of_fixed_from_second = 0; + std::vector> probing_values; +}; + template struct constraint_prop_t { constraint_prop_t(mip_solver_context_t& context); - bool apply_round( - solution_t& sol, - f_t lp_run_time_after_feasible, - timer_t& timer, - std::optional>> probing_candidates = std::nullopt); + bool apply_round(solution_t& sol, + f_t lp_run_time_after_feasible, + timer_t& timer, + std::optional>> + probing_config = std::nullopt); void sort_by_implied_slack_consumption(solution_t& sol, raft::device_span vars, bool problem_ii); @@ -52,14 +62,20 @@ struct constraint_prop_t { raft::device_span vars, std::mt19937 rng); - bool find_integer( - solution_t& sol, - solution_t& orig_sol, - f_t lp_run_time_after_feasible, - timer_t& timer, - std::optional>> probing_candidates = std::nullopt); + bool find_integer(solution_t& sol, + solution_t& orig_sol, + f_t lp_run_time_after_feasible, + timer_t& timer, + std::optional>> + probing_config = std::nullopt); void find_set_integer_vars(solution_t& sol, rmm::device_uvector& set_vars); void find_unset_integer_vars(solution_t& sol, rmm::device_uvector& set_vars); + thrust::pair generate_double_probing_pair( + const solution_t& sol, + const solution_t& orig_sol, + i_t unset_var_idx, + const std::optional>> probing_config, + bool bulk_rounding); std::tuple probing_values(const solution_t& sol, const solution_t& orig_sol, i_t idx); @@ -67,13 +83,10 @@ struct constraint_prop_t { void set_host_bounds(const solution_t& sol); bool probe(solution_t& sol, problem_t* original_problem, - const std::vector>& var_probe_val_pairs, + const std::tuple, std::vector, std::vector>& var_probe_vals, size_t* set_count_ptr, - rmm::device_uvector& unset_vars); - - void relax_crossing_bound_vars(solution_t& sol, - raft::device_span lower_bounds, - raft::device_span upper_bounds); + rmm::device_uvector& unset_vars, + std::optional>> probing_config); void collapse_crossing_bounds(problem_t& problem, problem_t& orig_problem, const raft::handle_t* handle_ptr); @@ -82,6 +95,11 @@ struct constraint_prop_t { void sort_by_frac(solution_t& sol, raft::device_span vars); void restore_bounds(solution_t& sol); void save_bounds(solution_t& sol); + void copy_bounds(rmm::device_uvector& output_lb, + rmm::device_uvector& output_ub, + const rmm::device_uvector& input_lb, + const rmm::device_uvector& input_ub, + const raft::handle_t* handle_ptr); void copy_bounds(rmm::device_uvector& output_lb, rmm::device_uvector& output_ub, rmm::device_uvector& output_assignment, @@ -90,11 +108,11 @@ struct constraint_prop_t { const rmm::device_uvector& input_assignment, const raft::handle_t* handle_ptr); void restore_original_bounds(solution_t& sol, solution_t& orig_sol); - std::vector> generate_bulk_rounding_vector( + std::tuple, std::vector, std::vector> generate_bulk_rounding_vector( const solution_t& sol, const solution_t& orig_sol, const std::vector& host_vars_to_set, - const std::optional>> probing_candidates); + const std::optional>> probing_config); bool is_problem_ii(problem_t& problem); void restore_original_bounds_on_unfixed(problem_t& problem, problem_t& original_problem, @@ -103,11 +121,17 @@ struct constraint_prop_t { problem_t& original_problem, timer_t& timer, const raft::handle_t* handle_ptr); - + bool handle_fixed_vars( + solution_t& sol, + problem_t* original_problem, + const std::tuple, std::vector, std::vector>& var_probe_vals, + size_t* set_count_ptr, + rmm::device_uvector& unset_vars); mip_solver_context_t& context; problem_t temp_problem; solution_t temp_sol; bound_presolve_t bounds_update; + multi_probe_t multi_probe; bounds_repair_t bounds_repair; rmm::device_uvector set_vars; rmm::device_uvector unset_vars; @@ -117,13 +141,16 @@ struct constraint_prop_t { rmm::device_uvector assignment_restore; std::vector curr_host_assignment; raft::random::PCGenerator rng; - bool recovery_mode = false; - bool rounding_ii = false; - i_t bounds_prop_interval = 1; - i_t n_iter_in_recovery = 0; + bool recovery_mode = false; + bool rounding_ii = false; + i_t selected_update = 0; + i_t bounds_prop_interval = 1; + i_t n_iter_in_recovery = 0; + i_t max_n_failed_repair_iterations = 1; timer_t max_timer{0.}; bool use_probing_cache = true; static repair_stats_t repair_stats; + bool single_rounding_only = false; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu b/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu index 93eb93693..b9408755c 100644 --- a/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu +++ b/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu @@ -951,13 +951,13 @@ bool lb_constraint_prop_t::find_integer( // if the constraint is not ii, run LP if (lb_bounds_update.infeas_constraints_count == 0 && !timeout_happened) { - run_lp_with_vars_fixed(*orig_sol.problem_ptr, - orig_sol, - orig_sol.problem_ptr->integer_indices, - context.settings.get_tolerances(), - context.lp_state, - lp_run_time_after_feasible, - true); + relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = lp_run_time_after_feasible; + lp_settings.tolerance = orig_sol.problem_ptr->tolerances.absolute_tolerance; + lp_settings.save_state = false; + lp_settings.return_first_feasible = true; + run_lp_with_vars_fixed( + *orig_sol.problem_ptr, orig_sol, orig_sol.problem_ptr->integer_indices, lp_settings); } bool res_feasible = orig_sol.compute_feasibility(); return res_feasible; diff --git a/cpp/src/mip/presolve/multi_probe.cu b/cpp/src/mip/presolve/multi_probe.cu index a40b5af6d..699a5f1dd 100644 --- a/cpp/src/mip/presolve/multi_probe.cu +++ b/cpp/src/mip/presolve/multi_probe.cu @@ -92,6 +92,8 @@ void multi_probe_t::resize(problem_t& problem) { upd_0.resize(problem); upd_1.resize(problem); + host_lb.resize(problem.n_variables); + host_ub.resize(problem.n_variables); } template @@ -315,6 +317,30 @@ termination_criterion_t multi_probe_t::bound_update_loop(problem_t +void multi_probe_t::update_device_bounds(const raft::handle_t* handle_ptr) +{ + cuopt_assert(upd_0.lb.size() == host_lb.size(), "size of variable lower bound mismatch"); + raft::copy(upd_0.lb.data(), host_lb.data(), upd_0.lb.size(), handle_ptr->get_stream()); + cuopt_assert(upd_0.ub.size() == host_ub.size(), "size of variable upper bound mismatch"); + raft::copy(upd_0.ub.data(), host_ub.data(), upd_0.ub.size(), handle_ptr->get_stream()); + cuopt_assert(upd_1.lb.size() == host_lb.size(), "size of variable lower bound mismatch"); + raft::copy(upd_1.lb.data(), host_lb.data(), upd_1.lb.size(), handle_ptr->get_stream()); + cuopt_assert(upd_1.ub.size() == host_ub.size(), "size of variable upper bound mismatch"); + raft::copy(upd_1.ub.data(), host_ub.data(), upd_1.ub.size(), handle_ptr->get_stream()); +} + +template +void multi_probe_t::update_host_bounds(const raft::handle_t* handle_ptr, + const raft::device_span variable_lb, + const raft::device_span variable_ub) +{ + cuopt_assert(variable_lb.size() == host_lb.size(), "size of variable lower bound mismatch"); + raft::copy(host_lb.data(), variable_lb.data(), variable_lb.size(), handle_ptr->get_stream()); + cuopt_assert(variable_ub.size() == host_ub.size(), "size of variable upper bound mismatch"); + raft::copy(host_ub.data(), variable_ub.data(), variable_ub.size(), handle_ptr->get_stream()); +} + template void multi_probe_t::copy_problem_into_probing_buffers(problem_t& pb, const raft::handle_t* handle_ptr) @@ -354,12 +380,16 @@ termination_criterion_t multi_probe_t::solve_for_interval( template termination_criterion_t multi_probe_t::solve( problem_t& pb, - const std::tuple, std::vector, std::vector>& var_probe_vals) + const std::tuple, std::vector, std::vector>& var_probe_vals, + bool use_host_bounds) { timer_t timer(settings.time_limit); auto& handle_ptr = pb.handle_ptr; - - copy_problem_into_probing_buffers(pb, handle_ptr); + if (use_host_bounds) { + update_device_bounds(handle_ptr); + } else { + copy_problem_into_probing_buffers(pb, handle_ptr); + } set_bounds(var_probe_vals, handle_ptr); return bound_update_loop(pb, handle_ptr, timer); diff --git a/cpp/src/mip/presolve/multi_probe.cuh b/cpp/src/mip/presolve/multi_probe.cuh index 8fb98c9b5..609702e90 100644 --- a/cpp/src/mip/presolve/multi_probe.cuh +++ b/cpp/src/mip/presolve/multi_probe.cuh @@ -42,7 +42,8 @@ class multi_probe_t { termination_criterion_t solve( problem_t& pb, - const std::tuple, std::vector, std::vector>& var_probe_vals); + const std::tuple, std::vector, std::vector>& var_probe_vals, + bool use_host_bounds = false); termination_criterion_t solve_for_interval( problem_t& pb, @@ -70,10 +71,15 @@ class multi_probe_t { const raft::handle_t* handle_ptr); void constraint_stats(problem_t& pb, const raft::handle_t* handle_ptr); void copy_problem_into_probing_buffers(problem_t& pb, const raft::handle_t* handle_ptr); - + void update_host_bounds(const raft::handle_t* handle_ptr, + const raft::device_span variable_lb, + const raft::device_span variable_ub); + void update_device_bounds(const raft::handle_t* handle_ptr); mip_solver_context_t& context; bounds_update_data_t upd_0; bounds_update_data_t upd_1; + std::vector host_lb; + std::vector host_ub; bool skip_0; bool skip_1; settings_t settings; diff --git a/cpp/src/mip/presolve/trivial_presolve.cuh b/cpp/src/mip/presolve/trivial_presolve.cuh index d7cfa6dc1..f4940a62f 100644 --- a/cpp/src/mip/presolve/trivial_presolve.cuh +++ b/cpp/src/mip/presolve/trivial_presolve.cuh @@ -277,9 +277,10 @@ void update_from_csr(problem_t& pb) pb.n_constraints = updated_n_cnst; pb.n_variables = updated_n_vars; - CUOPT_LOG_INFO("After trivial presolve updated number of %d constraints %d variables", + CUOPT_LOG_INFO("After trivial presolve updated %d constraints %d variables. Objective offset %f", updated_n_cnst, - updated_n_vars); + updated_n_vars, + pb.presolve_data.objective_offset); // check successive cnst in coo increases by atmost 1 // update csr offset pb.offsets.resize(pb.n_constraints + 1, handle_ptr->get_stream()); diff --git a/cpp/src/mip/problem/host_helper.cuh b/cpp/src/mip/problem/host_helper.cuh index fc3e3ee1b..4023a32d2 100644 --- a/cpp/src/mip/problem/host_helper.cuh +++ b/cpp/src/mip/problem/host_helper.cuh @@ -18,6 +18,9 @@ #pragma once #include + +#include + #include namespace cuopt::linear_programming::detail { diff --git a/cpp/src/mip/problem/problem.cu b/cpp/src/mip/problem/problem.cu index 494f5e40b..ea162cdfb 100644 --- a/cpp/src/mip/problem/problem.cu +++ b/cpp/src/mip/problem/problem.cu @@ -83,7 +83,6 @@ void problem_t::op_problem_cstr_body(const optimization_problem_t::problem_t( const typename mip_solver_settings_t::tolerances_t tolerances_) : original_problem_ptr(&problem_), handle_ptr(problem_.get_handle_ptr()), + integer_fixed_variable_map(problem_.get_n_variables(), problem_.get_handle_ptr()->get_stream()), + tolerances(tolerances_), n_variables(problem_.get_n_variables()), n_constraints(problem_.get_n_constraints()), n_binary_vars(0), @@ -139,7 +140,8 @@ problem_t::problem_t( var_names(problem_.get_variable_names()), row_names(problem_.get_row_names()), objective_name(problem_.get_objective_name()), - tolerances(tolerances_) + lp_state(*this, problem_.get_handle_ptr()->get_stream()), + fixing_helpers(n_constraints, n_variables, handle_ptr) { op_problem_cstr_body(problem_); branch_and_bound_callback = nullptr; @@ -150,6 +152,8 @@ problem_t::problem_t(const problem_t& problem_) : original_problem_ptr(problem_.original_problem_ptr), tolerances(problem_.tolerances), handle_ptr(problem_.handle_ptr), + integer_fixed_problem(problem_.integer_fixed_problem), + integer_fixed_variable_map(problem_.integer_fixed_variable_map, handle_ptr->get_stream()), branch_and_bound_callback(nullptr), n_variables(problem_.n_variables), n_constraints(problem_.n_constraints), @@ -185,7 +189,9 @@ problem_t::problem_t(const problem_t& problem_) row_names(problem_.row_names), objective_name(problem_.objective_name), is_scaled_(problem_.is_scaled_), - preprocess_called(problem_.preprocess_called) + preprocess_called(problem_.preprocess_called), + lp_state(problem_.lp_state), + fixing_helpers(problem_.fixing_helpers, handle_ptr) { } @@ -194,6 +200,8 @@ problem_t::problem_t(const problem_t& problem_, bool no_deep : original_problem_ptr(problem_.original_problem_ptr), tolerances(problem_.tolerances), handle_ptr(problem_.handle_ptr), + integer_fixed_problem(problem_.integer_fixed_problem), + integer_fixed_variable_map(problem_.n_variables, handle_ptr->get_stream()), n_variables(problem_.n_variables), n_constraints(problem_.n_constraints), n_binary_vars(problem_.n_binary_vars), @@ -279,7 +287,9 @@ problem_t::problem_t(const problem_t& problem_, bool no_deep row_names(problem_.row_names), objective_name(problem_.objective_name), is_scaled_(problem_.is_scaled_), - preprocess_called(problem_.preprocess_called) + preprocess_called(problem_.preprocess_called), + lp_state(problem_.lp_state), + fixing_helpers(problem_.fixing_helpers, handle_ptr) { } @@ -610,7 +620,6 @@ bool problem_t::pre_process_assignment(rmm::device_uvector& assig auto d_additional_var_id_per_var = cuopt::device_copy(presolve_data.additional_var_id_per_var, handle_ptr->get_stream()); - // handle free var logic by substituting the free vars and their corresponding vars thrust::for_each(handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), thrust::make_counting_iterator(original_problem_ptr->get_n_variables()), @@ -825,7 +834,7 @@ void problem_t::compute_related_variables(double time_limit) i_t slice_begin = i * max_slice_size; i_t slice_end = slice_begin + slice_size; - CUOPT_LOG_DEBUG("Iter %d: %d [%d %d] alloc'd %gmb", + CUOPT_LOG_TRACE("Iter %d: %d [%d %d] alloc'd %gmb", i, slice_size, slice_begin, @@ -1049,13 +1058,82 @@ void problem_t::fix_given_variables(problem_t& original_prob const rmm::device_uvector& variables_to_fix, const raft::handle_t* handle_ptr) { - i_t TPB = 64; - fix_given_variables_kernel<<get_stream()>>>( - original_problem.view(), - view(), - raft::device_span{assignment.data(), assignment.size()}, - raft::device_span{const_cast(variables_to_fix.data()), variables_to_fix.size()}); + fixing_helpers.reduction_in_rhs.resize(n_constraints, handle_ptr->get_stream()); + fixing_helpers.variable_fix_mask.resize(original_problem.n_variables, handle_ptr->get_stream()); + thrust::fill(handle_ptr->get_thrust_policy(), + fixing_helpers.reduction_in_rhs.begin(), + fixing_helpers.reduction_in_rhs.end(), + 0); + thrust::fill(handle_ptr->get_thrust_policy(), + fixing_helpers.variable_fix_mask.begin(), + fixing_helpers.variable_fix_mask.end(), + 0); + + thrust::for_each(handle_ptr->get_thrust_policy(), + variables_to_fix.begin(), + variables_to_fix.end(), + [variable_fix_mask = make_span(fixing_helpers.variable_fix_mask)] __device__( + i_t x) { variable_fix_mask[x] = 1; }); + const i_t num_segments = original_problem.n_constraints; + f_t initial_value{0.}; + + auto input_transform_it = thrust::make_transform_iterator( + thrust::make_counting_iterator(0), + [coefficients = make_span(original_problem.coefficients), + variables = make_span(original_problem.variables), + variable_fix_mask = make_span(fixing_helpers.variable_fix_mask), + assignment = make_span(assignment), + int_tol = original_problem.tolerances.integrality_tolerance] __device__(i_t idx) -> f_t { + i_t var_idx = variables[idx]; + if (variable_fix_mask[var_idx]) { + f_t reduction = coefficients[idx] * floor(assignment[var_idx] + int_tol); + return reduction; + } else { + return 0.; + } + }); + // Determine temporary device storage requirements + void* d_temp_storage = nullptr; + size_t temp_storage_bytes = 0; + cub::DeviceSegmentedReduce::Reduce(d_temp_storage, + temp_storage_bytes, + input_transform_it, + fixing_helpers.reduction_in_rhs.data(), + num_segments, + original_problem.offsets.data(), + original_problem.offsets.data() + 1, + cuda::std::plus<>{}, + initial_value, + handle_ptr->get_stream()); + + rmm::device_uvector temp_storage(temp_storage_bytes, handle_ptr->get_stream()); + d_temp_storage = thrust::raw_pointer_cast(temp_storage.data()); + + // Run reduction + cub::DeviceSegmentedReduce::Reduce(d_temp_storage, + temp_storage_bytes, + input_transform_it, + fixing_helpers.reduction_in_rhs.data(), + num_segments, + original_problem.offsets.data(), + original_problem.offsets.data() + 1, + cuda::std::plus<>{}, + initial_value, + handle_ptr->get_stream()); RAFT_CHECK_CUDA(handle_ptr->get_stream()); + thrust::for_each( + handle_ptr->get_thrust_policy(), + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(0) + n_constraints, + [lower_bounds = make_span(constraint_lower_bounds), + upper_bounds = make_span(constraint_upper_bounds), + original_lower_bounds = make_span(original_problem.constraint_lower_bounds), + original_upper_bounds = make_span(original_problem.constraint_upper_bounds), + reduction_in_rhs = make_span(fixing_helpers.reduction_in_rhs)] __device__(i_t cstr_idx) { + lower_bounds[cstr_idx] = original_lower_bounds[cstr_idx] - reduction_in_rhs[cstr_idx]; + upper_bounds[cstr_idx] = original_upper_bounds[cstr_idx] - reduction_in_rhs[cstr_idx]; + }); + handle_ptr->sync_stream(); } template @@ -1065,6 +1143,7 @@ problem_t problem_t::get_problem_after_fixing_vars( rmm::device_uvector& variable_map, const raft::handle_t* handle_ptr) { + auto start_time = std::chrono::high_resolution_clock::now(); cuopt_assert(n_variables == assignment.size(), "Assignment size issue"); problem_t problem(*this, true); CUOPT_LOG_DEBUG("Fixing %d variables", variables_to_fix.size()); @@ -1104,6 +1183,18 @@ problem_t problem_t::get_problem_after_fixing_vars( problem.reverse_original_ids[original_ids[h_variable_map[i]]] = i; } RAFT_CHECK_CUDA(handle_ptr->get_stream()); + auto end_time = std::chrono::high_resolution_clock::now(); + double time_taken = + std::chrono::duration_cast(end_time - start_time).count(); + static double total_time_taken = 0.; + static int total_calls = 0; + total_time_taken += time_taken; + total_calls++; + CUOPT_LOG_DEBUG( + "Time taken to fix variables: %f milliseconds, average: %f milliseconds total time: %f", + time_taken, + total_time_taken / total_calls, + total_time_taken); return problem; } @@ -1176,6 +1267,60 @@ void problem_t::remove_given_variables(problem_t& original_p check_problem_representation(true); } +template +rmm::device_uvector problem_t::get_fixed_assignment_from_integer_fixed_problem( + const rmm::device_uvector& assignment) +{ + rmm::device_uvector fixed_assignment(integer_fixed_variable_map.size(), + handle_ptr->get_stream()); + // first remove the assignment and variable related vectors + thrust::gather(handle_ptr->get_thrust_policy(), + integer_fixed_variable_map.begin(), + integer_fixed_variable_map.end(), + assignment.begin(), + fixed_assignment.begin()); + return fixed_assignment; +} + +template +void problem_t::compute_integer_fixed_problem() +{ + cuopt_assert(integer_fixed_problem == nullptr, "Integer fixed problem already computed"); + if (n_variables == n_integer_vars) { + integer_fixed_problem = nullptr; + return; + } + rmm::device_uvector assignment(n_variables, handle_ptr->get_stream()); + integer_fixed_problem = std::make_shared>(get_problem_after_fixing_vars( + assignment, integer_indices, integer_fixed_variable_map, handle_ptr)); + integer_fixed_problem->check_problem_representation(true); + integer_fixed_problem->lp_state.resize(*integer_fixed_problem, handle_ptr->get_stream()); +} + +template +void problem_t::copy_rhs_from_problem(const raft::handle_t* handle_ptr) +{ + raft::copy(integer_fixed_problem->constraint_lower_bounds.data(), + constraint_lower_bounds.data(), + integer_fixed_problem->constraint_lower_bounds.size(), + handle_ptr->get_stream()); + raft::copy(integer_fixed_problem->constraint_upper_bounds.data(), + constraint_upper_bounds.data(), + integer_fixed_problem->constraint_upper_bounds.size(), + handle_ptr->get_stream()); +} +template +void problem_t::fill_integer_fixed_problem(rmm::device_uvector& assignment, + const raft::handle_t* handle_ptr) +{ + cuopt_assert(integer_fixed_problem->n_variables > 0, "Integer fixed problem not computed"); + copy_rhs_from_problem(handle_ptr); + integer_fixed_problem->fix_given_variables(*this, assignment, integer_indices, handle_ptr); + combine_constraint_bounds(*integer_fixed_problem, + integer_fixed_problem->combined_bounds); + integer_fixed_problem->check_problem_representation(true); +} + template std::vector>> compute_var_to_constraint_map( const problem_t& pb) diff --git a/cpp/src/mip/problem/problem.cuh b/cpp/src/mip/problem/problem.cuh index b69168da7..cc29ce0e6 100644 --- a/cpp/src/mip/problem/problem.cuh +++ b/cpp/src/mip/problem/problem.cuh @@ -20,17 +20,20 @@ // THIS IS LIKELY THE INNER-MOST INCLUDE // FOR COMPILE TIME, WE SHOULD KEEP THE INCLUDES ON THIS HEADER MINIMAL +#include "host_helper.cuh" #include "presolve_data.cuh" +#include +#include + #include #include #include #include "host_helper.cuh" +#include "problem_fixing.cuh" #include -#include - #include #include #include @@ -96,7 +99,12 @@ class problem_t { void post_process_solution(solution_t& solution); void compute_transpose_of_problem(); f_t get_user_obj_from_solver_obj(f_t solver_obj); - + void compute_integer_fixed_problem(); + void fill_integer_fixed_problem(rmm::device_uvector& assignment, + const raft::handle_t* handle_ptr); + void copy_rhs_from_problem(const raft::handle_t* handle_ptr); + rmm::device_uvector get_fixed_assignment_from_integer_fixed_problem( + const rmm::device_uvector& assignment); bool is_integer(f_t val) const; bool integer_equal(f_t val1, f_t val2) const; @@ -195,6 +203,8 @@ class problem_t { const optimization_problem_t* original_problem_ptr; const raft::handle_t* handle_ptr; + std::shared_ptr> integer_fixed_problem = nullptr; + rmm::device_uvector integer_fixed_variable_map; std::function&)> branch_and_bound_callback; @@ -255,6 +265,13 @@ class problem_t { std::string objective_name; bool is_scaled_{false}; bool preprocess_called{false}; + // this LP state keeps the warm start data of some solution of + // 1. Original problem: it is unchanged and part of it is used + // to warm start slightly modified problems. + // 2. Integer fixed problem: this is useful as the problem structure + // is always the same and only the RHS changes. Using this helps in warm start. + lp_state_t lp_state; + problem_fixing_helpers_t fixing_helpers; }; } // namespace linear_programming::detail diff --git a/cpp/src/mip/problem/problem_fixing.cuh b/cpp/src/mip/problem/problem_fixing.cuh new file mode 100644 index 000000000..27c9dc84f --- /dev/null +++ b/cpp/src/mip/problem/problem_fixing.cuh @@ -0,0 +1,44 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#include + +namespace cuopt { +namespace linear_programming::detail { + +template +struct problem_fixing_helpers_t { + problem_fixing_helpers_t(i_t n_constraints, i_t n_variables, const raft::handle_t* handle_ptr) + : reduction_in_rhs(n_constraints, handle_ptr->get_stream()), + variable_fix_mask(n_variables, handle_ptr->get_stream()) + { + } + + problem_fixing_helpers_t(const problem_fixing_helpers_t& other, const raft::handle_t* handle_ptr) + : reduction_in_rhs(other.reduction_in_rhs, handle_ptr->get_stream()), + variable_fix_mask(other.variable_fix_mask, handle_ptr->get_stream()) + { + } + + rmm::device_uvector reduction_in_rhs; + rmm::device_uvector variable_fix_mask; +}; + +} // namespace linear_programming::detail +} // namespace cuopt \ No newline at end of file diff --git a/cpp/src/mip/problem/problem_helpers.cuh b/cpp/src/mip/problem/problem_helpers.cuh index d6bddfb13..dfcebc9d3 100644 --- a/cpp/src/mip/problem/problem_helpers.cuh +++ b/cpp/src/mip/problem/problem_helpers.cuh @@ -21,6 +21,7 @@ #include +#include #include #include @@ -179,6 +180,8 @@ __global__ void kernel_check_transpose_validity(raft::device_span coe __syncthreads(); // Would want to assert there but no easy way to gtest it, so moved it to the host if (!shared_found) { + DEVICE_LOG_DEBUG( + "For cstr %d, var %d, value %f was not found in the transpose", constraint_id, col, value); *failed = true; return; } diff --git a/cpp/src/mip/problem/problem_kernels.cuh b/cpp/src/mip/problem/problem_kernels.cuh index 71aeffe8f..bd25e4108 100644 --- a/cpp/src/mip/problem/problem_kernels.cuh +++ b/cpp/src/mip/problem/problem_kernels.cuh @@ -23,45 +23,6 @@ namespace cuopt::linear_programming::detail { -// TODO replace this kernel with actually removing the variables from the assignment -template -__global__ void fix_given_variables_kernel(const typename problem_t::view_t problem, - typename problem_t::view_t new_problem, - raft::device_span assignment, - raft::device_span vars_to_fix) -{ - // use a threadblock for each constraint - // loop over variables sequentially (assuming fixations occur with few variables) - i_t c = blockIdx.x; - const f_t int_tol = problem.tolerances.integrality_tolerance; - f_t th_total_reduction = 0.; - for (i_t v_idx = 0; v_idx < vars_to_fix.size(); ++v_idx) { - i_t var_to_fix = vars_to_fix[v_idx]; - f_t var_val = assignment[var_to_fix]; - auto [offset_begin, offset_end] = problem.range_for_constraint(c); - for (i_t i = threadIdx.x + offset_begin; i < offset_end; i += blockDim.x) { - i_t v = problem.variables[i]; - // search for var - if (v == var_to_fix) { - f_t coeff = problem.coefficients[i]; - f_t curr_val = floor(var_val + int_tol) * coeff; - // there should be only one var that is contained within the constraint - th_total_reduction += curr_val; - } - } - // wait until the var is fixed if any, so that at a time single var modifies the constraint - __syncthreads(); - } - __shared__ f_t shmem[raft::WarpSize]; - f_t total_reduction = raft::blockReduce(th_total_reduction, (char*)shmem); - // the new problem constraint values are not initialized, so we are basically initializng them - // here even if the total_reduction is zero, do the initiailizrion here - if (threadIdx.x == 0) { - new_problem.constraint_lower_bounds[c] = problem.constraint_lower_bounds[c] - total_reduction; - new_problem.constraint_upper_bounds[c] = problem.constraint_upper_bounds[c] - total_reduction; - } -} - template __global__ void compute_new_offsets(const typename problem_t::view_t orig_problem, typename problem_t::view_t new_problem, diff --git a/cpp/src/mip/relaxed_lp/lp_state.cuh b/cpp/src/mip/relaxed_lp/lp_state.cuh index e43c342be..0961a537f 100644 --- a/cpp/src/mip/relaxed_lp/lp_state.cuh +++ b/cpp/src/mip/relaxed_lp/lp_state.cuh @@ -22,6 +22,9 @@ namespace cuopt::linear_programming::detail { +template +class problem_t; + template class lp_state_t { public: diff --git a/cpp/src/mip/relaxed_lp/relaxed_lp.cu b/cpp/src/mip/relaxed_lp/relaxed_lp.cu index 05fcf3205..604f759df 100644 --- a/cpp/src/mip/relaxed_lp/relaxed_lp.cu +++ b/cpp/src/mip/relaxed_lp/relaxed_lp.cu @@ -33,20 +33,12 @@ namespace cuopt::linear_programming::detail { template -optimization_problem_solution_t get_relaxed_lp_solution(problem_t& op_problem, - solution_t& solution, - f_t tolerance, - f_t time_limit, - bool check_infeasibility, - bool return_first_feasible) +optimization_problem_solution_t get_relaxed_lp_solution( + problem_t& op_problem, + solution_t& solution, + const relaxed_lp_settings_t& settings) { - return get_relaxed_lp_solution(op_problem, - solution.assignment, - solution.lp_state, - tolerance, - time_limit, - check_infeasibility, - return_first_feasible); + return get_relaxed_lp_solution(op_problem, solution.assignment, solution.lp_state, settings); } template @@ -54,37 +46,28 @@ optimization_problem_solution_t get_relaxed_lp_solution( problem_t& op_problem, rmm::device_uvector& assignment, lp_state_t& lp_state, - f_t tolerance, - f_t time_limit, - bool check_infeasibility, - bool return_first_feasible, - bool save_state) + const relaxed_lp_settings_t& settings) { raft::common::nvtx::range fun_scope("get_relaxed_lp_solution"); - pdlp_solver_settings_t settings{}; - settings.detect_infeasibility = check_infeasibility; - settings.set_optimality_tolerance(tolerance); - settings.tolerances.relative_primal_tolerance = tolerance / 100.; - settings.tolerances.relative_dual_tolerance = tolerance / 100.; - settings.time_limit = time_limit; - if (return_first_feasible) { settings.per_constraint_residual = true; } - // settings.set_save_best_primal_so_far(true); - // currently disable first primal setting as it is not supported without per constraint mode - settings.first_primal_feasible = return_first_feasible; - pdlp_solver_t lp_solver(op_problem, settings); - if (save_state) { + pdlp_solver_settings_t pdlp_settings{}; + pdlp_settings.detect_infeasibility = settings.check_infeasibility; + pdlp_settings.set_optimality_tolerance(settings.tolerance); + pdlp_settings.tolerances.relative_primal_tolerance = settings.tolerance / 100.; + pdlp_settings.tolerances.relative_dual_tolerance = settings.tolerance / 100.; + pdlp_settings.time_limit = settings.time_limit; + if (settings.return_first_feasible) { pdlp_settings.per_constraint_residual = true; } + pdlp_settings.first_primal_feasible = settings.return_first_feasible; + pdlp_solver_t lp_solver(op_problem, pdlp_settings); + if (settings.save_state) { i_t prev_size = lp_state.prev_dual.size(); CUOPT_LOG_DEBUG( "setting initial primal solution of size %d dual size %d problem vars %d cstrs %d", - lp_state.prev_primal.size(), + assignment.size(), lp_state.prev_dual.size(), op_problem.n_variables, op_problem.n_constraints); lp_state.resize(op_problem, op_problem.handle_ptr->get_stream()); - raft::copy(lp_state.prev_primal.data(), - assignment.data(), - assignment.size(), - op_problem.handle_ptr->get_stream()); + clamp_within_var_bounds(assignment, &op_problem, op_problem.handle_ptr); // The previous dual sometimes contain invalid values w.r.t current problem // Adjust better dual values when we use warm start thrust::tabulate(op_problem.handle_ptr->get_thrust_policy(), @@ -95,7 +78,7 @@ optimization_problem_solution_t get_relaxed_lp_solution( if (!isfinite(x) || i >= prev_size) { return 0.0; } return x; }); - lp_solver.set_initial_primal_solution(lp_state.prev_primal); + lp_solver.set_initial_primal_solution(assignment); lp_solver.set_initial_dual_solution(lp_state.prev_dual); } CUOPT_LOG_DEBUG( @@ -108,7 +91,7 @@ optimization_problem_solution_t get_relaxed_lp_solution( auto solver_response = lp_solver.run_solver(start_time); if (solver_response.get_primal_solution().size() != 0 && - solver_response.get_dual_solution().size() != 0 && save_state) { + solver_response.get_dual_solution().size() != 0 && settings.save_state) { CUOPT_LOG_DEBUG("saving initial primal solution of size %d", lp_state.prev_primal.size()); lp_state.set_state(solver_response.get_primal_solution(), solver_response.get_dual_solution()); } @@ -129,24 +112,36 @@ optimization_problem_solution_t get_relaxed_lp_solution( return solver_response; } -// returns true if the problem is inevitablyinfeasible +// Run LP with variables fixed to specific values template bool run_lp_with_vars_fixed(problem_t& op_problem, + problem_t& fixed_problem, solution_t& solution, - const rmm::device_uvector& variables_to_fix, - typename mip_solver_settings_t::tolerances_t tols, - lp_state_t& lp_state, - f_t time_limit, - bool return_first_feasible, - bound_presolve_t* bound_presolve) + rmm::device_uvector& fixed_assignment, + rmm::device_uvector& variable_map, + relaxed_lp_settings_t& settings, + bound_presolve_t* bound_presolve, + bool check_fixed_assignment_feasibility) { - // if we are fixing all vars, there is no lp to be run - if (variables_to_fix.size() == (size_t)op_problem.n_variables) { return true; } - auto [fixed_problem, fixed_assignment, variable_map] = solution.fix_variables(variables_to_fix); + if (check_fixed_assignment_feasibility) { + solution_t temp_solution(fixed_problem); + raft::copy(temp_solution.assignment.data(), + fixed_assignment.data(), + fixed_assignment.size(), + fixed_problem.handle_ptr->get_stream()); + bool temp_solution_feasible = temp_solution.compute_feasibility(); + if (!temp_solution_feasible) { + CUOPT_LOG_DEBUG( + "Infeasible solution detected with fixed vars LP. Sol excess %f fixed sol excess %f", + solution.get_total_excess(), + temp_solution.get_total_excess()); + settings.time_limit = 1; + } + } if (bound_presolve != nullptr) { bound_presolve->resize(fixed_problem); // run bounds prop to quickly discover inevitably infeasible - bound_presolve->settings.time_limit = (time_limit / 10); + bound_presolve->settings.time_limit = (settings.time_limit / 10); auto term_crit = bound_presolve->solve(fixed_problem); bound_presolve->settings = {}; if (bound_presolve->infeas_constraints_count > 0) { @@ -157,50 +152,73 @@ bool run_lp_with_vars_fixed(problem_t& op_problem, } } fixed_problem.check_problem_representation(true); - const bool check_feas = false; // if we are on the original problem and fixing the integers, save the state // if we are in recombiners and on a smaller problem, don't update the state with integers fixed - bool save_state = false; - auto solver_response = get_relaxed_lp_solution(fixed_problem, - fixed_assignment, - lp_state, - tols.absolute_tolerance, - time_limit, - check_feas, - return_first_feasible, - save_state); + CUOPT_LOG_TRACE("save_state %d", settings.save_state); + auto& lp_state = fixed_problem.lp_state; + auto solver_response = + get_relaxed_lp_solution(fixed_problem, fixed_assignment, lp_state, settings); // unfix the assignment on given result no matter if it is feasible solution.unfix_variables(fixed_assignment, variable_map); if (bound_presolve != nullptr) { bound_presolve->resize(op_problem); } return false; } +// returns true if the problem is inevitably infeasible +template +bool run_lp_with_vars_fixed(problem_t& op_problem, + solution_t& solution, + const rmm::device_uvector& variables_to_fix, + relaxed_lp_settings_t& settings, + bound_presolve_t* bound_presolve, + bool check_fixed_assignment_feasibility, + bool use_integer_fixed_problem) +{ + // if we are fixing all vars, there is no lp to be run + if (variables_to_fix.size() == (size_t)op_problem.n_variables) { return true; } + if (use_integer_fixed_problem) { + op_problem.fill_integer_fixed_problem(solution.assignment, op_problem.handle_ptr); + auto fixed_assignment = + op_problem.get_fixed_assignment_from_integer_fixed_problem(solution.assignment); + return run_lp_with_vars_fixed(op_problem, + *op_problem.integer_fixed_problem, + solution, + fixed_assignment, + op_problem.integer_fixed_variable_map, + settings, + bound_presolve, + check_fixed_assignment_feasibility); + } else { + auto [fixed_problem, fixed_assignment, variable_map] = solution.fix_variables(variables_to_fix); + return run_lp_with_vars_fixed(op_problem, + fixed_problem, + solution, + fixed_assignment, + variable_map, + settings, + bound_presolve, + check_fixed_assignment_feasibility); + } +} + #define INSTANTIATE(F_TYPE) \ template optimization_problem_solution_t get_relaxed_lp_solution( \ problem_t & op_problem, \ solution_t & solution, \ - F_TYPE tolerance, \ - F_TYPE time_limit, \ - bool check_infeasibility, \ - bool return_first_feasible); \ + const relaxed_lp_settings_t& settings); \ template optimization_problem_solution_t get_relaxed_lp_solution( \ problem_t & op_problem, \ rmm::device_uvector & assignment, \ lp_state_t & lp_state, \ - F_TYPE tolerance, \ - F_TYPE time_limit, \ - bool check_infeasibility, \ - bool return_first_feasible, \ - bool save_state); \ + const relaxed_lp_settings_t& settings); \ template bool run_lp_with_vars_fixed( \ problem_t & op_problem, \ solution_t & solution, \ const rmm::device_uvector& variables_to_fix, \ - typename mip_solver_settings_t::tolerances_t tols, \ - lp_state_t& lp_state, \ - F_TYPE time_limit, \ - bool return_first_feasible, \ - bound_presolve_t* bound_presolve); + relaxed_lp_settings_t& settings, \ + bound_presolve_t* bound_presolve, \ + bool check_fixed_assignment_feasibility, \ + bool use_integer_fixed_problem); #if MIP_INSTANTIATE_FLOAT INSTANTIATE(float) diff --git a/cpp/src/mip/relaxed_lp/relaxed_lp.cuh b/cpp/src/mip/relaxed_lp/relaxed_lp.cuh index 39ccc7ef5..def1b7fa2 100644 --- a/cpp/src/mip/relaxed_lp/relaxed_lp.cuh +++ b/cpp/src/mip/relaxed_lp/relaxed_lp.cuh @@ -26,34 +26,35 @@ namespace cuopt::linear_programming::detail { +struct relaxed_lp_settings_t { + double tolerance = 1e-4; + double time_limit = 1.0; + bool check_infeasibility = true; + bool return_first_feasible = false; + bool save_state = true; + bool per_constraint_residual = false; +}; + template optimization_problem_solution_t get_relaxed_lp_solution( problem_t& op_problem, solution_t& solution, - f_t tolerance, - f_t time_limit = 20., - bool check_infeasibility = true, - bool return_first_feasible = false); + const relaxed_lp_settings_t& settings); template optimization_problem_solution_t get_relaxed_lp_solution( problem_t& op_problem, rmm::device_uvector& assignment, lp_state_t& lp_state, - f_t tolerance, - f_t time_limit = 20., - bool check_infeasibility = true, - bool return_first_feasible = false, - bool save_state = true); + const relaxed_lp_settings_t& settings); template bool run_lp_with_vars_fixed(problem_t& op_problem, solution_t& solution, const rmm::device_uvector& variables_to_fix, - typename mip_solver_settings_t::tolerances_t tols, - lp_state_t& lp_state, - f_t time_limit = 20., - bool return_first_feasible = false, - bound_presolve_t* bound_presolve = nullptr); + relaxed_lp_settings_t& settings, + bound_presolve_t* bound_presolve = nullptr, + bool check_fixed_assignment_feasibility = false, + bool use_integer_fixed_problem = false); } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/solution/solution.cu b/cpp/src/mip/solution/solution.cu index c4ca7be5b..74ec4c41c 100644 --- a/cpp/src/mip/solution/solution.cu +++ b/cpp/src/mip/solution/solution.cu @@ -310,7 +310,7 @@ bool solution_t::compute_feasibility() compute_infeasibility(); compute_number_of_integers(); i_t h_n_feas_constraints = n_feasible_constraints.value(handle_ptr->get_stream()); - is_feasible = h_n_feas_constraints == problem_ptr->n_constraints; + is_feasible = h_n_feas_constraints == problem_ptr->n_constraints && test_number_all_integer(); CUOPT_LOG_TRACE("is_feasible %d n_feasible_cstr %d all_cstr %d", is_feasible, h_n_feas_constraints, diff --git a/cpp/src/mip/solve.cu b/cpp/src/mip/solve.cu index 49ab7e5ba..dcfcdd0b1 100644 --- a/cpp/src/mip/solve.cu +++ b/cpp/src/mip/solve.cu @@ -126,7 +126,11 @@ mip_solution_t run_mip(detail::problem_t& problem, cuopt_func_call(auto saved_problem = scaled_problem); if (settings.mip_scaling) { scaling.scale_problem(); } - if (settings.has_initial_solution()) { scaling.scale_primal(settings.get_initial_solution()); } + if (settings.initial_solutions.size() > 0) { + for (const auto& initial_solution : settings.initial_solutions) { + scaling.scale_primal(*initial_solution); + } + } // only call preprocess on scaled problem, so we can compute feasibility on the original problem scaled_problem.preprocess_problem(); // cuopt_func_call((check_scaled_problem(scaled_problem, saved_problem))); diff --git a/cpp/src/mip/solver_context.cuh b/cpp/src/mip/solver_context.cuh index 3bb78380f..1f7f67623 100644 --- a/cpp/src/mip/solver_context.cuh +++ b/cpp/src/mip/solver_context.cuh @@ -33,11 +33,7 @@ struct mip_solver_context_t { problem_t* problem_ptr_, mip_solver_settings_t settings_, pdlp_initial_scaling_strategy_t& scaling) - : handle_ptr(handle_ptr_), - problem_ptr(problem_ptr_), - settings(settings_), - scaling(scaling), - lp_state(*problem_ptr) + : handle_ptr(handle_ptr_), problem_ptr(problem_ptr_), settings(settings_), scaling(scaling) { cuopt_assert(problem_ptr != nullptr, "problem_ptr is nullptr"); stats.solution_bound = problem_ptr->maximize ? std::numeric_limits::infinity() @@ -46,7 +42,6 @@ struct mip_solver_context_t { raft::handle_t const* const handle_ptr; problem_t* problem_ptr; - lp_state_t lp_state; const mip_solver_settings_t settings; pdlp_initial_scaling_strategy_t& scaling; solver_stats_t stats; diff --git a/cpp/src/mip/solver_settings.cu b/cpp/src/mip/solver_settings.cu index f6853776e..aae52caa4 100644 --- a/cpp/src/mip/solver_settings.cu +++ b/cpp/src/mip/solver_settings.cu @@ -23,17 +23,14 @@ namespace cuopt::linear_programming { template -void mip_solver_settings_t::set_initial_solution(const f_t* initial_solution, +void mip_solver_settings_t::add_initial_solution(const f_t* initial_solution, i_t size, rmm::cuda_stream_view stream) { cuopt_expects( initial_solution != nullptr, error_type_t::ValidationError, "initial_solution cannot be null"); - if (!initial_solution_) { - initial_solution_ = std::make_shared>(size, stream); - } - - raft::copy(initial_solution_.get()->data(), initial_solution, size, stream); + initial_solutions.emplace_back(std::make_shared>(size, stream)); + raft::copy(initial_solutions.back()->data(), initial_solution, size, stream); } template @@ -43,19 +40,6 @@ void mip_solver_settings_t::set_mip_callback( mip_callbacks_.push_back(callback); } -template -rmm::device_uvector& mip_solver_settings_t::get_initial_solution() const -{ - if (!initial_solution_) { throw std::runtime_error("Initial solution has not been set"); } - return *initial_solution_; -} - -template -bool mip_solver_settings_t::has_initial_solution() const -{ - return initial_solution_.get() != nullptr; -} - template const std::vector mip_solver_settings_t::get_mip_callbacks() const diff --git a/cpp/tests/mip/bounds_standardization_test.cu b/cpp/tests/mip/bounds_standardization_test.cu index 01b4a64ff..07e386bbd 100644 --- a/cpp/tests/mip/bounds_standardization_test.cu +++ b/cpp/tests/mip/bounds_standardization_test.cu @@ -76,10 +76,12 @@ void test_bounds_standardization_test(std::string test_instance) detail::solution_t solution_1(standardized_problem); mip_solver_settings_t default_settings{}; + detail::relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = 120.; + lp_settings.tolerance = default_settings.tolerances.absolute_tolerance; // run the problem through pdlp - auto result_1 = detail::get_relaxed_lp_solution( - standardized_problem, solution_1, default_settings.tolerances.absolute_tolerance, 120.); + auto result_1 = detail::get_relaxed_lp_solution(standardized_problem, solution_1, lp_settings); solution_1.compute_feasibility(); bool sol_1_feasible = (int)result_1.get_termination_status() == CUOPT_TERIMINATION_STATUS_OPTIMAL; // the problem might not be feasible in terms of per constraint residual diff --git a/cpp/tests/mip/elim_var_remap_test.cu b/cpp/tests/mip/elim_var_remap_test.cu index 4d1c2fd17..a1af26a33 100644 --- a/cpp/tests/mip/elim_var_remap_test.cu +++ b/cpp/tests/mip/elim_var_remap_test.cu @@ -161,9 +161,11 @@ void test_elim_var_solution(std::string test_instance) mip_solver_settings_t default_settings{}; detail::solution_t solution_1(standardized_problem); + detail::relaxed_lp_settings_t lp_settings; + lp_settings.time_limit = 120.; + lp_settings.tolerance = default_settings.tolerances.absolute_tolerance; // run the problem through pdlp - auto result_1 = detail::get_relaxed_lp_solution( - standardized_problem, solution_1, default_settings.tolerances.absolute_tolerance, 120.); + auto result_1 = detail::get_relaxed_lp_solution(standardized_problem, solution_1, lp_settings); solution_1.compute_feasibility(); // the solution might not be feasible per row as we are getting the result of pdlp bool sol_1_feasible = (int)result_1.get_termination_status() == CUOPT_TERIMINATION_STATUS_OPTIMAL; @@ -187,9 +189,11 @@ void test_elim_var_solution(std::string test_instance) trivial_presolve(sub_problem); detail::solution_t solution_2(sub_problem); + detail::relaxed_lp_settings_t lp_settings_2; + lp_settings_2.time_limit = 120.; + lp_settings_2.tolerance = default_settings.tolerances.absolute_tolerance; // run the problem through pdlp - auto result_2 = detail::get_relaxed_lp_solution( - sub_problem, solution_2, default_settings.tolerances.absolute_tolerance, 120.); + auto result_2 = detail::get_relaxed_lp_solution(sub_problem, solution_2, lp_settings_2); solution_2.compute_feasibility(); bool sol_2_feasible = (int)result_2.get_termination_status() == CUOPT_TERIMINATION_STATUS_OPTIMAL; EXPECT_EQ((int)result_2.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); diff --git a/python/cuopt/cuopt/linear_programming/solver/solver.pxd b/python/cuopt/cuopt/linear_programming/solver/solver.pxd index 3982e2ff8..255fcd764 100644 --- a/python/cuopt/cuopt/linear_programming/solver/solver.pxd +++ b/python/cuopt/cuopt/linear_programming/solver/solver.pxd @@ -84,7 +84,7 @@ cdef extern from "cuopt/linear_programming/solver_settings.hpp" namespace "cuopt ) except + # MIP settings - void set_initial_mip_solution( + void add_initial_mip_solution( const f_t* initial_solution, i_t size ) except + diff --git a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx index 9ba65bbc0..93a303489 100644 --- a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx +++ b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx @@ -305,7 +305,7 @@ cdef set_solver_setting( cdef uintptr_t callback_ptr = 0 if mip: if data_model_obj is not None and data_model_obj.get_initial_primal_solution().shape[0] != 0: # noqa - c_solver_settings.set_initial_mip_solution( + c_solver_settings.add_initial_mip_solution( c_initial_primal_solution, data_model_obj.get_initial_primal_solution().shape[0] )