diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 3163719488..713d55f16b 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -201,12 +201,12 @@ int run_single_file(std::string file_path, } } - settings.time_limit = time_limit; - settings.heuristics_only = heuristics_only; - settings.num_cpu_threads = num_cpu_threads; - settings.log_to_console = log_to_console; - // settings.tolerances.relative_tolerance = 1e-10; - // settings.tolerances.absolute_tolerance = 1e-6; + settings.time_limit = time_limit; + settings.heuristics_only = heuristics_only; + settings.num_cpu_threads = num_cpu_threads; + settings.log_to_console = log_to_console; + settings.tolerances.relative_tolerance = 1e-12; + settings.tolerances.absolute_tolerance = 1e-6; cuopt::linear_programming::benchmark_info_t benchmark_info; settings.benchmark_info_ptr = &benchmark_info; auto start_run_solver = std::chrono::high_resolution_clock::now(); diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 9067040129..ecd7488f19 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -422,6 +422,10 @@ branch_and_bound_t::branch_and_bound_t( template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { + // TODO remove this after submip PR merge + global_variables::mutex_upper.lock(); + global_variables::repair_queue.clear(); + global_variables::mutex_upper.unlock(); mip_status_t status = mip_status_t::UNSET; mip_solution_t incumbent(original_lp.num_cols); if (guess.size() != 0) { diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu index 5b9f60f44d..e0d96d447c 100644 --- a/cpp/src/mip/diversity/diversity_manager.cu +++ b/cpp/src/mip/diversity/diversity_manager.cu @@ -28,6 +28,7 @@ constexpr bool from_dir = false; constexpr bool fj_only_run = false; +constexpr bool fp_only_run = false; namespace cuopt::linear_programming::detail { @@ -112,8 +113,8 @@ bool diversity_manager_t::run_local_search(solution_t& solut timer_t& timer, ls_config_t& ls_config) { - // i_t ls_mab_option = mab_ls.select_mab_option(); - // mab_ls_config_t::get_local_search_and_lm_from_config(ls_mab_option, ls_config); + i_t ls_mab_option = mab_ls.select_mab_option(); + mab_ls_config_t::get_local_search_and_lm_from_config(ls_mab_option, ls_config); assignment_hash_map.insert(solution); constexpr i_t skip_solutions_threshold = 3; if (assignment_hash_map.check_skip_solution(solution, skip_solutions_threshold)) { return false; } @@ -293,7 +294,7 @@ void diversity_manager_t::generate_initial_solutions() population.var_threshold); population.print(); auto new_sol_vector = population.get_external_solutions(); - if (!fj_only_run) { recombine_and_ls_with_all(new_sol_vector); } + if (!fj_only_run && !fp_only_run) { recombine_and_ls_with_all(new_sol_vector); } } template @@ -322,7 +323,6 @@ bool diversity_manager_t::run_presolve(f_t time_limit) stats.presolve_time = presolve_timer.elapsed_time(); lp_optimal_solution.resize(problem_ptr->n_variables, problem_ptr->handle_ptr->get_stream()); problem_ptr->handle_ptr->sync_stream(); - cudaDeviceSynchronize(); return true; } @@ -382,6 +382,15 @@ void diversity_manager_t::run_fj_alone(solution_t& solution) CUOPT_LOG_INFO("FJ alone finished!"); } +// returns the best feasible solution +template +void diversity_manager_t::run_fp_alone(solution_t& solution) +{ + CUOPT_LOG_INFO("Running FP alone!"); + ls.run_fp(solution, timer, &population.weights, false); + CUOPT_LOG_INFO("FP alone finished!"); +} + // returns the best feasible solution template solution_t diversity_manager_t::run_solver() @@ -411,7 +420,7 @@ 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 - if (!from_dir) { generate_quick_feasible_solution(); } + if (!from_dir && !fp_only_run && !fj_only_run) { 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 = @@ -477,8 +486,10 @@ solution_t diversity_manager_t::run_solver() } population.allocate_solutions(); if (check_b_b_preemption()) { return population.best_feasible(); } - // generate a population with 5 solutions(FP+FJ) - generate_initial_solutions(); + if (!fp_only_run) { + // 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(); @@ -489,6 +500,12 @@ solution_t diversity_manager_t::run_solver() return population.best_feasible(); } + if (fp_only_run) { + auto sol = generate_solution(timer.remaining_time(), false); + run_fp_alone(sol); + return sol; + } + if (timer.check_time_limit()) { return population.best_feasible(); } main_loop(); return population.best_feasible(); diff --git a/cpp/src/mip/diversity/diversity_manager.cuh b/cpp/src/mip/diversity/diversity_manager.cuh index 712ad3f4a2..caf4d68036 100644 --- a/cpp/src/mip/diversity/diversity_manager.cuh +++ b/cpp/src/mip/diversity/diversity_manager.cuh @@ -46,6 +46,7 @@ class diversity_manager_t { // generates initial solutions void generate_initial_solutions(); void run_fj_alone(solution_t& solution); + void run_fp_alone(solution_t& solution); // main loop of diversity improvements void main_loop(); // randomly chooses a recombiner and returns the offspring diff --git a/cpp/src/mip/diversity/multi_armed_bandit.cuh b/cpp/src/mip/diversity/multi_armed_bandit.cuh index 52eef2c56b..d520f3202f 100644 --- a/cpp/src/mip/diversity/multi_armed_bandit.cuh +++ b/cpp/src/mip/diversity/multi_armed_bandit.cuh @@ -30,24 +30,13 @@ constexpr double ls_alpha = 0.03; template struct mab_ls_config_t { - static constexpr i_t n_of_ls = 2; - static constexpr i_t n_of_configs = 4; - static constexpr i_t n_of_arms = n_of_ls * n_of_configs; - static constexpr i_t ls_local_mins[n_of_configs] = {50, 100, 200, 500}; - static constexpr i_t ls_line_segment_local_mins[n_of_configs] = {10, 20, 40, 100}; + // RANDOM is the last ls method + static constexpr i_t n_of_arms = static_cast(ls_method_t::LS_METHODS_SIZE); static void get_local_search_and_lm_from_config(i_t config_id, ls_config_t& ls_config) { - ls_method_t local_search = ls_method_t(config_id % n_of_ls); - i_t lm_id = config_id / n_of_ls; - if (local_search == ls_method_t::FJ_LINE_SEGMENT) { - ls_config.ls_method = ls_method_t::FJ_LINE_SEGMENT; - ls_config.n_local_mins_for_line_segment = ls_line_segment_local_mins[lm_id]; - } else { - ls_config.ls_method = ls_method_t::FJ_ANNEALING; - ls_config.n_local_mins = ls_local_mins[lm_id]; - } - mab_ls_config_t::last_lm_config = lm_id; + ls_method_t local_search = ls_method_t(config_id % n_of_arms); + ls_config.ls_method = local_search; mab_ls_config_t::last_ls_mab_option = config_id; } static i_t last_lm_config; @@ -61,11 +50,9 @@ i_t mab_ls_config_t::last_ls_mab_option = 0; struct ls_work_normalized_reward_t { int option_id; - static constexpr double reward_per_option[mab_ls_config_t::n_of_configs] = { - 2, 1, 0.5, 0.25}; ls_work_normalized_reward_t(int option_id) : option_id(option_id) {} - double operator()(double factor) const { return factor * reward_per_option[option_id]; } + double operator()(double factor) const { return factor; } }; struct recombiner_work_normalized_reward_t { diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump.cu b/cpp/src/mip/feasibility_jump/feasibility_jump.cu index fa88bddd19..8e864dcf20 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump.cu +++ b/cpp/src/mip/feasibility_jump/feasibility_jump.cu @@ -239,13 +239,31 @@ void fj_t::populate_climber_views() template void fj_t::copy_weights(const weight_t& weights, - const raft::handle_t* handle_ptr) + const raft::handle_t* handle_ptr, + std::optional new_size) { - cuopt_assert(cstr_weights.size() == weights.cstr_weights.size(), "Size mismatch"); - raft::copy(cstr_weights.data(), - weights.cstr_weights.data(), - weights.cstr_weights.size(), - handle_ptr->get_stream()); + i_t old_size = weights.cstr_weights.size(); + cstr_weights.resize(new_size.value_or(weights.cstr_weights.size()), handle_ptr->get_stream()); + cstr_left_weights.resize(new_size.value_or(weights.cstr_weights.size()), + handle_ptr->get_stream()); + cstr_right_weights.resize(new_size.value_or(weights.cstr_weights.size()), + handle_ptr->get_stream()); + thrust::for_each(handle_ptr->get_thrust_policy(), + thrust::counting_iterator(0), + thrust::counting_iterator(new_size.value_or(weights.cstr_weights.size())), + [old_size, + fj_weights = make_span(cstr_weights), + fj_left_weights = make_span(cstr_left_weights), + fj_right_weights = make_span(cstr_right_weights), + new_weights = make_span(weights.cstr_weights)] __device__(i_t idx) { + fj_weights[idx] = idx >= old_size ? 1. : new_weights[idx]; + // TODO: ask Alice how we can manage the previous left,right weights + fj_left_weights[idx] = idx >= old_size ? 1. : new_weights[idx]; + fj_right_weights[idx] = idx >= old_size ? 1. : new_weights[idx]; + cuopt_assert(isfinite(fj_weights[idx]), "invalid weight"); + cuopt_assert(isfinite(fj_left_weights[idx]), "invalid left weight"); + cuopt_assert(isfinite(fj_right_weights[idx]), "invalid right weight"); + }); thrust::transform(handle_ptr->get_thrust_policy(), weights.objective_weight.data(), weights.objective_weight.data() + 1, diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump.cuh b/cpp/src/mip/feasibility_jump/feasibility_jump.cuh index 028b1c8baf..5d362a3d3d 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump.cuh +++ b/cpp/src/mip/feasibility_jump/feasibility_jump.cuh @@ -205,7 +205,9 @@ class fj_t { void set_fj_settings(fj_settings_t settings_); void reset_weights(const rmm::cuda_stream_view& stream, f_t weight = 10.); void randomize_weights(const raft::handle_t* handle_ptr); - void copy_weights(const weight_t& weights, const raft::handle_t* handle_ptr); + void copy_weights(const weight_t& weights, + const raft::handle_t* handle_ptr, + std::optional new_size = std::nullopt); i_t host_loop(solution_t& solution, i_t climber_idx = 0); void run_step_device(i_t climber_idx = 0, bool use_graph = true); void run_step_device(const rmm::cuda_stream_view& stream, diff --git a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu index ee04155d6c..8286d8148e 100644 --- a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu +++ b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu @@ -288,7 +288,7 @@ bool feasibility_pump_t::run_fj_cycle_escape(solution_t& sol bool is_feasible; fj.settings.mode = fj_mode_t::EXIT_NON_IMPROVING; fj.settings.update_weights = true; - fj.settings.feasibility_run = true; + fj.settings.feasibility_run = false; fj.settings.n_of_minimums_for_exit = 5000; fj.settings.time_limit = std::min(3., timer.remaining_time()); is_feasible = fj.solve(solution); @@ -344,9 +344,10 @@ bool feasibility_pump_t::handle_cycle(solution_t& solution) cuopt_assert(solution.test_number_all_integer(), "All must be integers before fj"); is_feasible = run_fj_cycle_escape(solution); cuopt_assert(solution.test_number_all_integer(), "All must be integers after fj"); - if (!is_feasible && cycle_queue.check_cycle(solution)) { + if (cycle_queue.check_cycle(solution)) { CUOPT_LOG_DEBUG("Cycle couldn't be broken. Perturbating FP"); perturbate(solution); + is_feasible = solution.get_feasible(); } cycle_queue.n_iterations_without_cycle = 0; cycle_queue.update_recent_solutions(solution); @@ -381,7 +382,6 @@ void feasibility_pump_t::reset() total_fp_time_until_cycle = 0; fp_fj_cycle_time_begin = timer.remaining_time(); max_n_of_integers = 0; - run_intensive_restart = false; config.alpha = default_alpha; last_distances.resize(0); } diff --git a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cuh b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cuh index 8b7891b04b..60a249f893 100644 --- a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cuh +++ b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cuh @@ -169,8 +169,7 @@ class feasibility_pump_t { f_t proj_and_round_time; f_t proj_begin; i_t n_fj_single_descents; - i_t max_n_of_integers = 0; - bool run_intensive_restart = false; + i_t max_n_of_integers = 0; cuopt::timer_t timer; }; diff --git a/cpp/src/mip/local_search/local_search.cu b/cpp/src/mip/local_search/local_search.cu index 711612550b..5e39491264 100644 --- a/cpp/src/mip/local_search/local_search.cu +++ b/cpp/src/mip/local_search/local_search.cu @@ -49,7 +49,8 @@ local_search_t::local_search_t(mip_solver_context_t& context lb_constraint_prop, line_segment_search, lp_optimal_solution_), - rng(cuopt::seed_generator::get_seed()) + rng(cuopt::seed_generator::get_seed()), + problem_with_objective_cut(*context.problem_ptr, context.problem_ptr->handle_ptr) { } @@ -105,17 +106,24 @@ bool local_search_t::run_local_search(solution_t& solution, fj_settings.update_weights = false; fj_settings.feasibility_run = false; fj.set_fj_settings(fj_settings); - fj.copy_weights(weights, solution.handle_ptr); - bool is_feas; - i_t rd = std::uniform_int_distribution(0, 1)(rng); + bool is_feas = false; + ls_method_t rd = static_cast(std::uniform_int_distribution( + static_cast(ls_method_t::FJ_LINE_SEGMENT), static_cast(ls_method_t::FP_SEARCH))(rng)); if (ls_config.ls_method == ls_method_t::FJ_LINE_SEGMENT) { rd = ls_method_t::FJ_LINE_SEGMENT; } else if (ls_config.ls_method == ls_method_t::FJ_ANNEALING) { rd = ls_method_t::FJ_ANNEALING; + } else if (ls_config.ls_method == ls_method_t::FP_SEARCH) { + rd = ls_method_t::FP_SEARCH; } if (rd == ls_method_t::FJ_LINE_SEGMENT && lp_optimal_exists) { + fj.copy_weights(weights, solution.handle_ptr); is_feas = run_fj_line_segment(solution, timer, ls_config); + } else if (rd == ls_method_t::FP_SEARCH) { + timer = timer_t(std::min(3., timer.remaining_time())); + is_feas = run_fp(solution, timer, &weights, false); } else { + fj.copy_weights(weights, solution.handle_ptr); is_feas = run_fj_annealing(solution, timer, ls_config); } return is_feas; @@ -140,7 +148,6 @@ bool local_search_t::run_fj_until_timer(solution_t& solution return is_feasible; } -// SIMULATED ANNEALING not fully implemented yet, placeholder template bool local_search_t::run_fj_annealing(solution_t& solution, timer_t timer, @@ -190,10 +197,12 @@ bool local_search_t::check_fj_on_lp_optimal(solution_t& solu bool perturb, timer_t timer) { - raft::copy(solution.assignment.data(), - lp_optimal_solution.data(), - solution.assignment.size(), - solution.handle_ptr->get_stream()); + if (lp_optimal_exists) { + raft::copy(solution.assignment.data(), + lp_optimal_solution.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + } if (perturb) { CUOPT_LOG_DEBUG("Perturbating solution on initial fj on optimal run!"); f_t perturbation_ratio = 0.2; @@ -218,7 +227,7 @@ bool local_search_t::check_fj_on_lp_optimal(solution_t& solu fj.settings.mode = fj_mode_t::EXIT_NON_IMPROVING; fj.settings.n_of_minimums_for_exit = 20000; fj.settings.update_weights = true; - fj.settings.feasibility_run = true; + fj.settings.feasibility_run = false; fj.settings.time_limit = std::min(30., timer.remaining_time()); fj.solve(solution); return solution.get_feasible(); @@ -235,7 +244,7 @@ bool local_search_t::run_fj_on_zero(solution_t& solution, ti fj.settings.mode = fj_mode_t::EXIT_NON_IMPROVING; fj.settings.n_of_minimums_for_exit = 20000; fj.settings.update_weights = true; - fj.settings.feasibility_run = true; + fj.settings.feasibility_run = false; fj.settings.time_limit = std::min(30., timer.remaining_time()); bool is_feasible = fj.solve(solution); return is_feasible; @@ -319,32 +328,124 @@ void local_search_t::resize_vectors(problem_t& problem, const raft::handle_t* handle_ptr) { fj_sol_on_lp_opt.resize(problem.n_variables, handle_ptr->get_stream()); - fj.resize_vectors(handle_ptr); - // TODO resize LNS when it is used fp.resize_vectors(problem, handle_ptr); } template -bool local_search_t::run_fp(solution_t& solution, timer_t timer) +void save_best_fp_solution(solution_t& solution, + rmm::device_uvector& best_solution, + f_t& best_objective, + bool feasibility_run) +{ + if (feasibility_run || solution.get_objective() < best_objective) { + CUOPT_LOG_DEBUG("Found better feasible in FP with obj %f. Continue with FJ!", + solution.get_objective()); + best_objective = solution.get_objective(); + raft::copy(best_solution.data(), + solution.assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + solution.problem_ptr->add_cutting_plane_at_objective(solution.get_objective() - + OBJECTIVE_EPSILON); + } +} + +template +void local_search_t::save_solution_and_add_cutting_plane( + solution_t& solution, rmm::device_uvector& best_solution, f_t& best_objective) +{ + if (solution.get_objective() < best_objective) { + raft::copy(best_solution.data(), + solution.assignment.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + best_objective = solution.get_objective(); + solution.problem_ptr->add_cutting_plane_at_objective(solution.get_objective() - + OBJECTIVE_EPSILON); + } +} + +template +bool local_search_t::run_fp(solution_t& solution, + timer_t timer, + const weight_t* weights, + bool feasibility_run) { const i_t n_fp_iterations = 1000000; - bool is_feasible = false; + bool is_feasible = solution.compute_feasibility(); + double best_objective = solution.get_objective(); + rmm::device_uvector best_solution(solution.assignment, solution.handle_ptr->get_stream()); + problem_t* old_problem_ptr = solution.problem_ptr; + fp.timer = timer_t(timer.remaining_time()); + if (!feasibility_run) { + // if it has not been initialized yet, create a new problem and move it to the cut problem + if (!problem_with_objective_cut.cutting_plane_added) { + problem_with_objective_cut = std::move(problem_t(*old_problem_ptr)); + } + problem_with_objective_cut.add_cutting_plane_at_objective(solution.get_objective() - + OBJECTIVE_EPSILON); + // Do the copy here for proper handling of the added constraints weight + fj.copy_weights(*weights, solution.handle_ptr, problem_with_objective_cut.n_constraints); + solution.problem_ptr = &problem_with_objective_cut; + solution.resize_to_problem(); + resize_vectors(problem_with_objective_cut, solution.handle_ptr); + lb_constraint_prop.temp_problem.setup(problem_with_objective_cut); + lb_constraint_prop.bounds_update.setup(lb_constraint_prop.temp_problem); + constraint_prop.bounds_update.resize(problem_with_objective_cut); + } for (i_t i = 0; i < n_fp_iterations && !timer.check_time_limit(); ++i) { - if (timer.check_time_limit()) { return false; } + if (timer.check_time_limit()) { + is_feasible = false; + break; + } CUOPT_LOG_DEBUG("fp_loop it %d", i); is_feasible = fp.run_single_fp_descent(solution); // if feasible return true if (is_feasible) { - return true; + if (feasibility_run) { + is_feasible = true; + break; + } else { + CUOPT_LOG_DEBUG("Found feasible in FP with obj %f. Continue with FJ!", + solution.get_objective()); + save_solution_and_add_cutting_plane(solution, best_solution, best_objective); + fp.config.alpha = default_alpha; + } } // if not feasible, it means it is a cycle else { - if (timer.check_time_limit()) { return false; } + if (timer.check_time_limit()) { + is_feasible = false; + break; + } is_feasible = fp.restart_fp(solution); - if (is_feasible) { return true; } + if (is_feasible) { + if (feasibility_run) { + is_feasible = true; + break; + } else { + CUOPT_LOG_DEBUG("Found feasible in FP with obj %f. Continue with FJ!", + solution.get_objective()); + save_solution_and_add_cutting_plane(solution, best_solution, best_objective); + fp.config.alpha = default_alpha; + } + } } } - return false; + if (!feasibility_run) { + raft::copy(solution.assignment.data(), + best_solution.data(), + solution.assignment.size(), + solution.handle_ptr->get_stream()); + solution.problem_ptr = old_problem_ptr; + solution.resize_to_problem(); + resize_vectors(*old_problem_ptr, solution.handle_ptr); + lb_constraint_prop.temp_problem.setup(*old_problem_ptr); + lb_constraint_prop.bounds_update.setup(lb_constraint_prop.temp_problem); + constraint_prop.bounds_update.resize(*old_problem_ptr); + solution.handle_ptr->sync_stream(); + } + return is_feasible; } template diff --git a/cpp/src/mip/local_search/local_search.cuh b/cpp/src/mip/local_search/local_search.cuh index fdb3ff5803..a2664e890b 100644 --- a/cpp/src/mip/local_search/local_search.cuh +++ b/cpp/src/mip/local_search/local_search.cuh @@ -25,7 +25,14 @@ namespace cuopt::linear_programming::detail { -enum ls_method_t { FJ_ANNEALING = 0, FJ_LINE_SEGMENT, RANDOM }; +// make sure RANDOM is always the last +enum class ls_method_t : int { + FJ_ANNEALING = 0, + FJ_LINE_SEGMENT, + FP_SEARCH, + RANDOM, + LS_METHODS_SIZE = RANDOM +}; template struct ls_config_t { @@ -66,8 +73,14 @@ class local_search_t { 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); bool run_staged_fp(solution_t& solution, timer_t timer, bool& early_exit); - bool run_fp(solution_t& solution, timer_t timer); + bool run_fp(solution_t& solution, + timer_t timer, + const weight_t* weights = nullptr, + bool feasibility_run = true); void resize_vectors(problem_t& problem, const raft::handle_t* handle_ptr); + void save_solution_and_add_cutting_plane(solution_t& solution, + rmm::device_uvector& best_solution, + f_t& best_objective); mip_solver_context_t& context; rmm::device_uvector& lp_optimal_solution; @@ -80,6 +93,7 @@ class local_search_t { line_segment_search_t line_segment_search; feasibility_pump_t fp; std::mt19937 rng; + problem_t problem_with_objective_cut; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/problem/problem.cu b/cpp/src/mip/problem/problem.cu index 888388ef18..b3ca958d2d 100644 --- a/cpp/src/mip/problem/problem.cu +++ b/cpp/src/mip/problem/problem.cu @@ -92,8 +92,8 @@ void problem_t::op_problem_cstr_body(const optimization_problem_tget_stream()); compute_n_integer_vars(); compute_binary_var_table(); + compute_vars_with_objective_coeffs(); } - compute_transpose_of_problem(); // Check after modifications check_problem_representation(true, is_mip); @@ -193,7 +193,8 @@ problem_t::problem_t(const problem_t& problem_) is_scaled_(problem_.is_scaled_), preprocess_called(problem_.preprocess_called), lp_state(problem_.lp_state), - fixing_helpers(problem_.fixing_helpers, handle_ptr) + fixing_helpers(problem_.fixing_helpers, handle_ptr), + vars_with_objective_coeffs(problem_.vars_with_objective_coeffs) { } @@ -291,7 +292,8 @@ problem_t::problem_t(const problem_t& problem_, bool no_deep is_scaled_(problem_.is_scaled_), preprocess_called(problem_.preprocess_called), lp_state(problem_.lp_state), - fixing_helpers(problem_.fixing_helpers, handle_ptr) + fixing_helpers(problem_.fixing_helpers, handle_ptr), + vars_with_objective_coeffs(problem_.vars_with_objective_coeffs) { } @@ -716,7 +718,7 @@ void problem_t::recompute_auxilliary_data(bool check_representation) { compute_n_integer_vars(); compute_binary_var_table(); - + compute_vars_with_objective_coeffs(); // TODO: speedup compute related variables const double time_limit = 30.; compute_related_variables(time_limit); @@ -999,6 +1001,7 @@ void problem_t::insert_variables(variables_delta_t& h_vars) compute_n_integer_vars(); compute_binary_var_table(); + compute_vars_with_objective_coeffs(); } // note that these don't change the reverse structure @@ -1590,6 +1593,42 @@ f_t problem_t::get_user_obj_from_solver_obj(f_t solver_obj) return presolve_data.objective_scaling_factor * (solver_obj + presolve_data.objective_offset); } +template +void problem_t::compute_vars_with_objective_coeffs() +{ + auto h_objective_coefficients = cuopt::host_copy(objective_coefficients); + std::vector vars_with_objective_coeffs_; + std::vector objective_coeffs_; + for (i_t i = 0; i < n_variables; ++i) { + if (h_objective_coefficients[i] != 0) { + vars_with_objective_coeffs_.push_back(i); + objective_coeffs_.push_back(h_objective_coefficients[i]); + } + } + vars_with_objective_coeffs = std::make_pair(vars_with_objective_coeffs_, objective_coeffs_); +} + +template +void problem_t::add_cutting_plane_at_objective(f_t objective) +{ + CUOPT_LOG_INFO("Adding cutting plane at objective %f", objective); + if (cutting_plane_added) { + // modify the RHS + i_t last_constraint = n_constraints - 1; + constraint_upper_bounds.set_element_async(last_constraint, objective, handle_ptr->get_stream()); + return; + } + cutting_plane_added = true; + constraints_delta_t h_constraints; + h_constraints.add_constraint(vars_with_objective_coeffs.first, + vars_with_objective_coeffs.second, + -std::numeric_limits::infinity(), + objective); + insert_constraints(h_constraints); + compute_transpose_of_problem(); + cuopt_func_call(check_problem_representation(true)); +} + #if MIP_INSTANTIATE_FLOAT template class problem_t; #endif diff --git a/cpp/src/mip/problem/problem.cuh b/cpp/src/mip/problem/problem.cuh index cc29ce0e65..d8b57ac6f3 100644 --- a/cpp/src/mip/problem/problem.cuh +++ b/cpp/src/mip/problem/problem.cuh @@ -112,6 +112,8 @@ class problem_t { cuopt::linear_programming::dual_simplex::user_problem_t& user_problem) const; void write_as_mps(const std::string& path); + void add_cutting_plane_at_objective(f_t objective); + void compute_vars_with_objective_coeffs(); struct view_t { DI std::pair reverse_range_for_var(i_t v) const @@ -272,6 +274,8 @@ class problem_t { // 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; + bool cutting_plane_added{false}; + std::pair, std::vector> vars_with_objective_coeffs; }; } // namespace linear_programming::detail diff --git a/cpp/src/mip/solve.cu b/cpp/src/mip/solve.cu index 7071c33fea..33cb40f70a 100644 --- a/cpp/src/mip/solve.cu +++ b/cpp/src/mip/solve.cu @@ -155,8 +155,8 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, mip_solver_settings_t const& settings) { try { - const f_t time_limit = - settings.time_limit == 0 ? std::numeric_limits::max() : settings.time_limit; + constexpr f_t max_time_limit = 1000000000; + const f_t time_limit = settings.time_limit == 0 ? max_time_limit : settings.time_limit; if (settings.heuristics_only && time_limit == std::numeric_limits::max()) { CUOPT_LOG_ERROR("Time limit cannot be infinity when heuristics only is set"); cuopt_expects(false,