diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 5880b3f924..ed9728cc94 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -103,7 +103,8 @@ bool check_guess(const lp_problem_t& original_lp, const std::vector& guess, f_t& primal_error, f_t& bound_error, - i_t& num_fractional) + i_t& num_fractional, + f_t& integrality_error) { bool feasible = false; std::vector residual(original_lp.num_rows); @@ -111,8 +112,13 @@ bool check_guess(const lp_problem_t& original_lp, matrix_vector_multiply(original_lp.A, 1.0, guess, -1.0, residual); primal_error = vector_norm_inf(residual); bound_error = 0.0; + integrality_error = 0.0; constexpr bool verbose = false; for (i_t j = 0; j < original_lp.num_cols; j++) { + if (var_types[j] == variable_type_t::INTEGER) { + integrality_error = std::max(integrality_error, std::abs(guess[j] - std::round(guess[j]))); + } + // l_j <= x_j infeas means x_j < l_j or l_j - x_j > 0 const f_t low_bound_err = std::max(0.0, original_lp.lower[j] - guess[j]); // x_j <= u_j infeas means u_j < x_j or x_j - u_j > 0 @@ -460,6 +466,64 @@ void branch_and_bound_t::update_user_bound(const lp_problem_t +void branch_and_bound_t::snap_to_lattice(const lp_problem_t& lp, + mip_node_t* node_ptr, + f_t leaf_obj) +{ + if (lp.objective_step.has_step()) { + f_t step = lp.objective_step.step_size; + f_t bias = lp.objective_step.bias; + // Round up to next value on the lattice: k * step + bias >= leaf_obj + f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); + node_ptr->lower_bound = k * step + bias; + } else if (lp.objective_is_integral) { + node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); + } +} + +template +f_t branch_and_bound_t::compute_node_cutoff(const lp_problem_t& lp, + f_t upper_bound) const +{ + if (lp.objective_step.has_step()) { + f_t step = lp.objective_step.step_size; + f_t bias = lp.objective_step.bias; + f_t k = std::ceil((upper_bound - bias) / step - settings_.integer_tol); + return (k - 1) * step + bias + settings_.dual_tol; + } else if (lp.objective_is_integral) { + return std::ceil(upper_bound - settings_.integer_tol) - 1 + settings_.dual_tol; + } + return upper_bound + settings_.dual_tol; +} + +template +bool branch_and_bound_t::accept_new_incumbent(f_t obj) +{ + if (!incumbent_.has_incumbent) return true; + + f_t threshold = incumbent_.objective - settings_.primal_tol; + f_t step = 0.0; + if (original_lp_.objective_step.has_step()) { + step = original_lp_.objective_step.step_size; + } else if (original_lp_.objective_is_integral) { + step = 1.0; + } + if (step > 0.0) { + // Half a step above the greatest lattice value strictly below the incumbent. + const f_t lattice_below = + compute_node_cutoff(original_lp_, incumbent_.objective) - settings_.dual_tol; + threshold = std::min(threshold, lattice_below + 0.5 * step); + } + return obj < threshold; +} + template bool branch_and_bound_t::set_solution_from_heuristics(const std::vector& solution, heuristics_origin_t origin) @@ -483,9 +547,10 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto feasible_solution_symbol(origin), compute_user_objective(original_lp_, obj)); - if (!incumbent_.has_incumbent || obj < incumbent_.objective) { + if (accept_new_incumbent(obj)) { f_t primal_err; f_t bound_err; + f_t integral_err; i_t num_fractional; mutex_original_lp_.lock(); if (crushed_solution.size() != original_lp_.num_cols) { @@ -494,19 +559,27 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto crush_primal_solution( original_problem_, original_lp_, solution, new_slacks_, crushed_solution); } - is_feasible = check_guess( - original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); + is_feasible = check_guess(original_lp_, + settings_, + var_types_, + crushed_solution, + primal_err, + bound_err, + num_fractional, + integral_err); mutex_original_lp_.unlock(); - mutex_upper_.lock(); - if (is_feasible && (!incumbent_.has_incumbent || obj < incumbent_.objective)) { - f_t current_upper_bound = upper_bound_.load(); - upper_bound_ = std::min(current_upper_bound, obj); - incumbent_.set_incumbent_solution(obj, crushed_solution); - if (current_upper_bound > upper_bound_.load()) { - report_heuristic(obj, origin); - success = true; - } + if (is_feasible) { + std::lock_guard lock(mutex_upper_); + if (accept_new_incumbent(obj)) { + f_t current_upper_bound = upper_bound_.load(); + upper_bound_ = std::min(current_upper_bound, obj); + incumbent_.set_incumbent_solution(obj, crushed_solution); + if (current_upper_bound > upper_bound_.load()) { + report_heuristic(obj, origin); + success = true; + } + } } else { attempt_repair = true; constexpr bool verbose = false; @@ -519,7 +592,6 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto num_fractional); } } - mutex_upper_.unlock(); } else { settings_.log.debug("Solution objective not better than current upper_bound_. Not accepted.\n"); } @@ -569,8 +641,15 @@ void branch_and_bound_t::queue_external_solution_deterministic( f_t primal_err; f_t bound_err; i_t num_fractional; - bool is_feasible = check_guess( - original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); + f_t integral_err; + bool is_feasible = check_guess(original_lp_, + settings_, + var_types_, + crushed_solution, + primal_err, + bound_err, + num_fractional, + integral_err); mutex_original_lp_.unlock(); if (!is_feasible) { @@ -697,13 +776,15 @@ bool branch_and_bound_t::repair_solution(const std::vector& edge_ f_t primal_error; f_t bound_error; i_t num_fractional; + f_t integral_err; feasible = check_guess(original_lp_, settings_, var_types_, lp_solution.x, primal_error, bound_error, - num_fractional); + num_fractional, + integral_err); repaired_obj = compute_objective(original_lp_, repaired_solution); constexpr bool verbose = false; if (verbose) { @@ -745,8 +826,7 @@ void branch_and_bound_t::repair_heuristic_solutions() repair_solution(edge_norms_, crushed_solution, repaired_obj, repaired_solution); if (is_feasible) { mutex_upper_.lock(); - - if (!incumbent_.has_incumbent || repaired_obj < incumbent_.objective) { + if (accept_new_incumbent(repaired_obj)) { upper_bound_ = std::min(upper_bound_.load(), repaired_obj); incumbent_.set_incumbent_solution(repaired_obj, repaired_solution); report_heuristic(repaired_obj, heuristics_origin_t::HEURISTICS); @@ -920,7 +1000,7 @@ void branch_and_bound_t::add_feasible_solution(const lp_problem_t -void branch_and_bound_t::snap_to_lattice(const lp_problem_t& lp, - mip_node_t* node_ptr, - f_t leaf_obj) -{ - if (lp.objective_step.has_step()) { - f_t step = lp.objective_step.step_size; - f_t bias = lp.objective_step.bias; - // Round up to next value on the lattice: k * step + bias >= leaf_obj - f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); - node_ptr->lower_bound = k * step + bias; - } else if (lp.objective_is_integral) { - node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); - } -} - template template std::pair branch_and_bound_t::update_tree_impl( @@ -1581,21 +1639,7 @@ dual_status_t branch_and_bound_t::solve_node_lp( lp_settings.concurrent_halt = &node_concurrent_halt_; } lp_settings.set_log(false); - f_t cutoff = upper_bound_.load(); - if (worker->leaf_problem.objective_step.has_step()) { - f_t step = worker->leaf_problem.objective_step.step_size; - f_t bias = worker->leaf_problem.objective_step.bias; - // Any improving feasible solution must have objective <= cutoff - step. - f_t k = std::floor((cutoff - bias) / step + settings_.integer_tol); - lp_settings.cut_off = (k - 1) * step + bias + settings_.dual_tol; - } else if (worker->leaf_problem.objective_is_integral) { - // If the objective is integral, any feasible solution should produce an upper bound that is - // (approximately) integral. We add a small tolerance and floor this value to get an integer, - // we then subtract 1, to stop simplex on problems that cannot improve the primal objective. - lp_settings.cut_off = std::floor(cutoff + settings_.integer_tol) - 1 + settings_.dual_tol; - } else { - lp_settings.cut_off = cutoff + settings_.dual_tol; - } + lp_settings.cut_off = compute_node_cutoff(worker->leaf_problem, upper_bound_.load()); lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); if (lp_settings.time_limit <= 0.0) { return dual_status_t::TIME_LIMIT; } @@ -4033,8 +4077,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t primal_err; f_t bound_err; i_t num_fractional; - const bool feasible = check_guess( - original_lp_, settings_, var_types_, crushed_guess, primal_err, bound_err, num_fractional); + f_t integral_err; + const bool feasible = check_guess(original_lp_, + settings_, + var_types_, + crushed_guess, + primal_err, + bound_err, + num_fractional, + integral_err); if (feasible) { const f_t computed_obj = compute_objective(original_lp_, crushed_guess); mutex_upper_.lock(); @@ -5378,7 +5429,7 @@ void branch_and_bound_t::deterministic_process_worker_solutions( deterministic_current_horizon_); bool improved = false; - if (!incumbent_.has_incumbent || sol->objective < incumbent_.objective) { + if (accept_new_incumbent(sol->objective)) { upper_bound_ = std::min(upper_bound_.load(), sol->objective); incumbent_.set_incumbent_solution(sol->objective, sol->solution); current_upper = sol->objective; @@ -5541,8 +5592,7 @@ void branch_and_bound_t::deterministic_sort_replay_events( // Process heuristic solution at its correct work unit timestamp position f_t new_upper = std::numeric_limits::infinity(); - - if (!incumbent_.has_incumbent || hsol.objective < incumbent_.objective) { + if (accept_new_incumbent(hsol.objective)) { upper_bound_ = std::min(upper_bound_.load(), hsol.objective); incumbent_.set_incumbent_solution(hsol.objective, hsol.solution); new_upper = hsol.objective; diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index db62fa745c..c626c44356 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -360,6 +360,9 @@ class branch_and_bound_t { mip_node_t* node_ptr, f_t leaf_obj); + f_t compute_node_cutoff(const simplex::lp_problem_t& lp, f_t upper_bound) const; + bool accept_new_incumbent(f_t obj); + // Launch a new best-first worker from a given bfs worker. void launch_bfs_worker(bfs_worker_t* worker); diff --git a/cpp/src/mip_heuristics/solution/solution.cu b/cpp/src/mip_heuristics/solution/solution.cu index 0d411a0ffc..6198b0e4c5 100644 --- a/cpp/src/mip_heuristics/solution/solution.cu +++ b/cpp/src/mip_heuristics/solution/solution.cu @@ -567,13 +567,10 @@ f_t solution_t::compute_max_int_violation() { return thrust::transform_reduce( handle_ptr->get_thrust_policy(), - thrust::make_counting_iterator(0), - thrust::make_counting_iterator(0) + problem_ptr->n_integer_vars, + problem_ptr->integer_indices.begin(), + problem_ptr->integer_indices.end(), cuda::proclaim_return_type([v = view()] __device__(i_t idx) -> f_t { - if (v.problem.variable_types[idx] == var_t::INTEGER) { - return abs(v.assignment[idx] - round(v.assignment[idx])); - } - return 0.; + return abs(v.assignment[idx] - round(v.assignment[idx])); }), 0., thrust::maximum()); diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 9558beda9d..e5613f2484 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -847,9 +847,19 @@ mip_solution_t solve_mip_helper( full_sol.copy_new_assignment( cuopt::host_copy(primal_solution, op_problem.get_handle_ptr()->get_stream())); full_sol.compute_feasibility(); + CUOPT_LOG_DEBUG( + "Post solve: max residual error=%g, max integrality error=%g, max variable bound " + "error=%g", + full_sol.compute_max_constraint_violation(), + full_sol.compute_max_int_violation(), + full_sol.compute_max_variable_violation()); if (!full_sol.get_feasible()) { CUOPT_LOG_WARN("The solution is not feasible after post solve"); } + if (mip::has_variable_bounds_violation( + full_problem.handle_ptr, full_sol.assignment, &full_problem)) { + CUOPT_LOG_WARN("The solution violates the variable bounds after post solve"); + } auto full_stats = sol.get_stats(); // add third party presolve time to cuopt presolve time