diff --git a/cpp/src/linear_programming/initial_scaling_strategy/initial_scaling.cu b/cpp/src/linear_programming/initial_scaling_strategy/initial_scaling.cu index 4c6cbf475b..72931267ad 100644 --- a/cpp/src/linear_programming/initial_scaling_strategy/initial_scaling.cu +++ b/cpp/src/linear_programming/initial_scaling_strategy/initial_scaling.cu @@ -375,18 +375,13 @@ void pdlp_initial_scaling_strategy_t::scale_problem() primal_size_h_, stream_view_); - raft::linalg::eltwiseDivideCheckZero( - const_cast&>(op_problem_scaled_.variable_lower_bounds).data(), - op_problem_scaled_.variable_lower_bounds.data(), - cummulative_variable_scaling_.data(), - primal_size_h_, - stream_view_); - raft::linalg::eltwiseDivideCheckZero( - const_cast&>(op_problem_scaled_.variable_upper_bounds).data(), - op_problem_scaled_.variable_upper_bounds.data(), - cummulative_variable_scaling_.data(), - primal_size_h_, - stream_view_); + using f_t2 = typename type_2::type; + cub::DeviceTransform::Transform(cuda::std::make_tuple(op_problem_scaled_.variable_bounds.data(), + cummulative_variable_scaling_.data()), + op_problem_scaled_.variable_bounds.data(), + primal_size_h_, + divide_check_zero(), + stream_view_); raft::linalg::eltwiseMultiply( const_cast&>(op_problem_scaled_.constraint_lower_bounds).data(), diff --git a/cpp/src/linear_programming/pdhg.cu b/cpp/src/linear_programming/pdhg.cu index f932eeb8d8..38373391aa 100644 --- a/cpp/src/linear_programming/pdhg.cu +++ b/cpp/src/linear_programming/pdhg.cu @@ -143,18 +143,18 @@ void pdhg_solver_t::compute_primal_projection_with_gradient( // project by max(min(x[i], upperbound[i]),lowerbound[i]) // compute delta_primal x'-x + using f_t2 = typename type_2::type; // All is fused in a single call to limit number of read / write in memory cub::DeviceTransform::Transform( cuda::std::make_tuple(current_saddle_point_state_.get_primal_solution().data(), problem_ptr->objective_coefficients.data(), current_saddle_point_state_.get_current_AtY().data(), - problem_ptr->variable_lower_bounds.data(), - problem_ptr->variable_upper_bounds.data()), + problem_ptr->variable_bounds.data()), thrust::make_zip_iterator(potential_next_primal_solution_.data(), current_saddle_point_state_.get_delta_primal().data(), tmp_primal_.data()), primal_size_h_, - primal_projection(primal_step_size.data()), + primal_projection(primal_step_size.data()), stream_view_); } diff --git a/cpp/src/linear_programming/pdlp.cu b/cpp/src/linear_programming/pdlp.cu index 84b04d43f9..824539e0ca 100644 --- a/cpp/src/linear_programming/pdlp.cu +++ b/cpp/src/linear_programming/pdlp.cu @@ -1041,20 +1041,21 @@ optimization_problem_solution_t pdlp_solver_t::run_solver( // Project initial primal solution if (pdlp_hyper_params::project_initial_primal) { - raft::linalg::ternaryOp(pdhg_solver_.get_primal_solution().data(), - pdhg_solver_.get_primal_solution().data(), - op_problem_scaled_.variable_lower_bounds.data(), - op_problem_scaled_.variable_upper_bounds.data(), - primal_size_h_, - clamp(), - stream_view_); - raft::linalg::ternaryOp(unscaled_primal_avg_solution_.data(), - unscaled_primal_avg_solution_.data(), - op_problem_scaled_.variable_lower_bounds.data(), - op_problem_scaled_.variable_upper_bounds.data(), - primal_size_h_, - clamp(), - stream_view_); + using f_t2 = typename type_2::type; + cub::DeviceTransform::Transform( + cuda::std::make_tuple(pdhg_solver_.get_primal_solution().data(), + op_problem_scaled_.variable_bounds.data()), + pdhg_solver_.get_primal_solution().data(), + primal_size_h_, + clamp(), + stream_view_); + cub::DeviceTransform::Transform( + cuda::std::make_tuple(unscaled_primal_avg_solution_.data(), + op_problem_scaled_.variable_bounds.data()), + unscaled_primal_avg_solution_.data(), + primal_size_h_, + clamp(), + stream_view_); } if (verbose) { diff --git a/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu b/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu index 236ec373d8..cf23c8a1d7 100644 --- a/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu +++ b/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu @@ -1388,6 +1388,14 @@ median breakpoint is evaluated, eliminating half of the components. The process is iterated until the argmin is identified. */ +template +struct extract_bounds_t { + __device__ thrust::tuple operator()(f_t2 bounds) + { + return thrust::make_tuple(get_lower(bounds), get_upper(bounds)); + } +}; + template void pdlp_restart_strategy_t::solve_bound_constrained_trust_region( localized_duality_gap_container_t& duality_gap, @@ -1451,10 +1459,13 @@ void pdlp_restart_strategy_t::solve_bound_constrained_trust_region( // component becomes fixed by its bounds // Copying primal / dual bound before sorting them according to threshold - raft::copy( - lower_bound_.data(), problem_ptr->variable_lower_bounds.data(), primal_size_h_, stream_view_); - raft::copy( - upper_bound_.data(), problem_ptr->variable_upper_bounds.data(), primal_size_h_, stream_view_); + using f_t2 = typename type_2::type; + cub::DeviceTransform::Transform( + problem_ptr->variable_bounds.data(), + thrust::make_zip_iterator(thrust::make_tuple(lower_bound_.data(), upper_bound_.data())), + primal_size_h_, + extract_bounds_t(), + stream_view_); raft::copy(lower_bound_.data() + primal_size_h_, transformed_constraint_lower_bounds_.data(), dual_size_h_, @@ -1632,13 +1643,13 @@ void pdlp_restart_strategy_t::solve_bound_constrained_trust_region( a_add_scalar_times_b(target_threshold_.data()), stream_view_); // project by max(min(x[i], upperbound[i]),lowerbound[i]) for primal part - raft::linalg::ternaryOp(duality_gap.primal_solution_tr_.data(), - duality_gap.primal_solution_tr_.data(), - problem_ptr->variable_lower_bounds.data(), - problem_ptr->variable_upper_bounds.data(), - primal_size_h_, - clamp(), - stream_view_); + using f_t2 = typename type_2::type; + cub::DeviceTransform::Transform(cuda::std::make_tuple(duality_gap.primal_solution_tr_.data(), + problem_ptr->variable_bounds.data()), + duality_gap.primal_solution_tr_.data(), + primal_size_h_, + clamp(), + stream_view_); // project by max(min(y[i], upperbound[i]),lowerbound[i]) raft::linalg::ternaryOp(duality_gap.dual_solution_tr_.data(), @@ -1646,7 +1657,7 @@ void pdlp_restart_strategy_t::solve_bound_constrained_trust_region( transformed_constraint_lower_bounds_.data(), transformed_constraint_upper_bounds_.data(), dual_size_h_, - clamp(), + constraint_clamp(), stream_view_); // } } diff --git a/cpp/src/linear_programming/termination_strategy/convergence_information.cu b/cpp/src/linear_programming/termination_strategy/convergence_information.cu index 2378b8b9bf..ce579e3104 100644 --- a/cpp/src/linear_programming/termination_strategy/convergence_information.cu +++ b/cpp/src/linear_programming/termination_strategy/convergence_information.cu @@ -336,7 +336,7 @@ void convergence_information_t::compute_dual_objective( problem_ptr->constraint_lower_bounds.data(), problem_ptr->constraint_upper_bounds.data(), dual_size_h_, - bound_value_reduced_cost_product(), + constraint_bound_value_reduced_cost_product(), stream_view_); cub::DeviceReduce::Sum(rmm_tmp_buffer_.data(), @@ -371,13 +371,13 @@ void convergence_information_t::compute_reduced_cost_from_primal_gradi { raft::common::nvtx::range fun_scope("compute_reduced_cost_from_primal_gradient"); - raft::linalg::ternaryOp(bound_value_.data(), - primal_gradient.data(), - problem_ptr->variable_lower_bounds.data(), - problem_ptr->variable_upper_bounds.data(), - primal_size_h_, - bound_value_gradient(), - stream_view_); + using f_t2 = typename type_2::type; + cub::DeviceTransform::Transform( + cuda::std::make_tuple(primal_gradient.data(), problem_ptr->variable_bounds.data()), + bound_value_.data(), + primal_size_h_, + bound_value_gradient(), + stream_view_); if (pdlp_hyper_params::handle_some_primal_gradients_on_finite_bounds_as_residuals) { raft::linalg::ternaryOp(reduced_cost_.data(), @@ -402,15 +402,15 @@ void convergence_information_t::compute_reduced_costs_dual_objective_c { raft::common::nvtx::range fun_scope("compute_reduced_costs_dual_objective_contribution"); + using f_t2 = typename type_2::type; // if reduced cost is positive -> lower bound, negative -> upper bounds, 0 -> 0 // if bound_val is not finite let element be -inf, otherwise bound_value*reduced_cost - raft::linalg::ternaryOp(bound_value_.data(), - reduced_cost_.data(), - problem_ptr->variable_lower_bounds.data(), - problem_ptr->variable_upper_bounds.data(), - primal_size_h_, - bound_value_reduced_cost_product(), - stream_view_); + cub::DeviceTransform::Transform( + cuda::std::make_tuple(reduced_cost_.data(), problem_ptr->variable_bounds.data()), + bound_value_.data(), + primal_size_h_, + bound_value_reduced_cost_product(), + stream_view_); // sum over bound_value*reduced_cost, but should be -inf if any element is -inf cub::DeviceReduce::Sum(rmm_tmp_buffer_.data(), diff --git a/cpp/src/linear_programming/termination_strategy/infeasibility_information.cu b/cpp/src/linear_programming/termination_strategy/infeasibility_information.cu index fd9651d9a2..1e1bf3abf2 100644 --- a/cpp/src/linear_programming/termination_strategy/infeasibility_information.cu +++ b/cpp/src/linear_programming/termination_strategy/infeasibility_information.cu @@ -254,16 +254,14 @@ void infeasibility_information_t::compute_max_violation( // Convert raw pointer to thrust::device_ptr to write directly device side through reduce thrust::device_ptr primal_ray_max_violation(primal_ray_max_violation_.data()); + using f_t2 = typename type_2::type; *primal_ray_max_violation = thrust::transform_reduce( handle_ptr_->get_thrust_policy(), - thrust::make_zip_iterator(thrust::make_tuple(primal_ray.data(), - problem_ptr->variable_lower_bounds.data(), - problem_ptr->variable_upper_bounds.data())), thrust::make_zip_iterator( - thrust::make_tuple(primal_ray.data() + primal_size_h_, - problem_ptr->variable_lower_bounds.data() + primal_size_h_, - problem_ptr->variable_upper_bounds.data() + primal_size_h_)), - max_violation(), + thrust::make_tuple(primal_ray.data(), problem_ptr->variable_bounds.data())), + thrust::make_zip_iterator(thrust::make_tuple( + primal_ray.data() + primal_size_h_, problem_ptr->variable_bounds.data() + primal_size_h_)), + max_violation(), f_t(0.0), thrust::maximum()); } @@ -329,7 +327,7 @@ void infeasibility_information_t::compute_homogenous_dual_objective( problem_ptr->constraint_lower_bounds.data(), problem_ptr->constraint_upper_bounds.data(), dual_size_h_, - bound_value_reduced_cost_product(), + constraint_bound_value_reduced_cost_product(), stream_view_); cub::DeviceReduce::Sum(rmm_tmp_buffer_.data(), @@ -364,13 +362,13 @@ template void infeasibility_information_t::compute_reduced_cost_from_primal_gradient( rmm::device_uvector& primal_gradient, rmm::device_uvector& primal_ray) { - raft::linalg::ternaryOp(bound_value_.data(), - primal_gradient.data(), - problem_ptr->variable_lower_bounds.data(), - problem_ptr->variable_upper_bounds.data(), - primal_size_h_, - bound_value_gradient(), - stream_view_); + using f_t2 = typename type_2::type; + cub::DeviceTransform::Transform( + cuda::std::make_tuple(primal_gradient.data(), problem_ptr->variable_bounds.data()), + bound_value_.data(), + primal_size_h_, + bound_value_gradient(), + stream_view_); if (pdlp_hyper_params::handle_some_primal_gradients_on_finite_bounds_as_residuals) { raft::linalg::ternaryOp(reduced_cost_.data(), @@ -393,16 +391,16 @@ void infeasibility_information_t::compute_reduced_cost_from_primal_gra template void infeasibility_information_t::compute_reduced_costs_dual_objective_contribution() { + using f_t2 = typename type_2::type; // Check if these bounds are the same as computed above // if reduced cost is positive -> lower bound, negative -> upper bounds, 0 -> 0 // if bound_val is not finite let element be -inf, otherwise bound_value*reduced_cost - raft::linalg::ternaryOp(bound_value_.data(), - reduced_cost_.data(), - problem_ptr->variable_lower_bounds.data(), - problem_ptr->variable_upper_bounds.data(), - primal_size_h_, - bound_value_reduced_cost_product(), - stream_view_); + cub::DeviceTransform::Transform( + cuda::std::make_tuple(reduced_cost_.data(), problem_ptr->variable_bounds.data()), + bound_value_.data(), + primal_size_h_, + bound_value_reduced_cost_product(), + stream_view_); // sum over bound_value*reduced_cost cub::DeviceReduce::Sum(rmm_tmp_buffer_.data(), diff --git a/cpp/src/linear_programming/translate.hpp b/cpp/src/linear_programming/translate.hpp index 6731be6ec3..0f1fc9f1c1 100644 --- a/cpp/src/linear_programming/translate.hpp +++ b/cpp/src/linear_programming/translate.hpp @@ -81,9 +81,9 @@ static dual_simplex::user_problem_t cuopt_problem_to_simplex_problem( } } user_problem.num_range_rows = user_problem.range_rows.size(); - user_problem.lower = cuopt::host_copy(model.variable_lower_bounds); - user_problem.upper = cuopt::host_copy(model.variable_upper_bounds); - user_problem.problem_name = model.original_problem_ptr->get_problem_name(); + std::tie(user_problem.lower, user_problem.upper) = + extract_host_bounds(model.variable_bounds, model.handle_ptr); + user_problem.problem_name = model.original_problem_ptr->get_problem_name(); if (model.row_names.size() > 0) { user_problem.row_names.resize(m); for (int i = 0; i < m; ++i) { @@ -164,8 +164,7 @@ void translate_to_crossover_problem(const detail::problem_t& problem, lp.obj_constant = problem.presolve_data.objective_offset; lp.obj_scale = problem.presolve_data.objective_scaling_factor; - std::vector lower = cuopt::host_copy(problem.variable_lower_bounds); - std::vector upper = cuopt::host_copy(problem.variable_upper_bounds); + auto [lower, upper] = extract_host_bounds(problem.variable_bounds, problem.handle_ptr); std::vector constraint_lower = cuopt::host_copy(problem.constraint_lower_bounds); std::vector constraint_upper = cuopt::host_copy(problem.constraint_upper_bounds); diff --git a/cpp/src/linear_programming/utilities/problem_checking.cu b/cpp/src/linear_programming/utilities/problem_checking.cu index d0fc6811b3..b53021a959 100644 --- a/cpp/src/linear_programming/utilities/problem_checking.cu +++ b/cpp/src/linear_programming/utilities/problem_checking.cu @@ -263,21 +263,16 @@ template void problem_checking_t::check_scaled_problem( detail::problem_t const& scaled_problem, detail::problem_t const& op_problem) { + using f_t2 = typename type_2::type; // original problem to host - auto& d_variable_upper_bounds = op_problem.variable_upper_bounds; - auto& d_variable_lower_bounds = op_problem.variable_lower_bounds; - auto& d_variable_types = op_problem.variable_types; - std::vector variable_upper_bounds(d_variable_upper_bounds.size()); - std::vector variable_lower_bounds(d_variable_lower_bounds.size()); + auto& d_variable_bounds = op_problem.variable_bounds; + auto& d_variable_types = op_problem.variable_types; + std::vector variable_bounds(d_variable_bounds.size()); std::vector variable_types(d_variable_types.size()); - raft::copy(variable_upper_bounds.data(), - d_variable_upper_bounds.data(), - d_variable_upper_bounds.size(), - op_problem.handle_ptr->get_stream()); - raft::copy(variable_lower_bounds.data(), - d_variable_lower_bounds.data(), - d_variable_lower_bounds.size(), + raft::copy(variable_bounds.data(), + d_variable_bounds.data(), + d_variable_bounds.size(), op_problem.handle_ptr->get_stream()); raft::copy(variable_types.data(), d_variable_types.data(), @@ -285,26 +280,20 @@ void problem_checking_t::check_scaled_problem( op_problem.handle_ptr->get_stream()); // scaled problem to host - std::vector scaled_variable_upper_bounds(scaled_problem.variable_upper_bounds.size()); - std::vector scaled_variable_lower_bounds(scaled_problem.variable_lower_bounds.size()); - std::vector scaled_variables(scaled_problem.variable_lower_bounds.size()); + std::vector scaled_variable_bounds(scaled_problem.variable_bounds.size()); - raft::copy(scaled_variable_upper_bounds.data(), - scaled_problem.variable_upper_bounds.data(), - scaled_problem.variable_upper_bounds.size(), - op_problem.handle_ptr->get_stream()); - raft::copy(scaled_variable_lower_bounds.data(), - scaled_problem.variable_lower_bounds.data(), - scaled_problem.variable_lower_bounds.size(), + raft::copy(scaled_variable_bounds.data(), + scaled_problem.variable_bounds.data(), + scaled_problem.variable_bounds.size(), op_problem.handle_ptr->get_stream()); for (size_t i = 0; i < variable_types.size(); ++i) { auto var_type = variable_types[i]; if (var_type == var_t::INTEGER) { // Integers should be untouched - cuopt_assert(variable_upper_bounds[i] == scaled_variable_upper_bounds[i], - "Mismatch upper scaling"); - cuopt_assert(variable_lower_bounds[i] == scaled_variable_lower_bounds[i], + cuopt_assert(get_lower(variable_bounds[i]) == get_lower(scaled_variable_bounds[i]), "Mismatch lower scaling"); + cuopt_assert(get_upper(variable_bounds[i]) == get_upper(scaled_variable_bounds[i]), + "Mismatch upper scaling"); } } } @@ -313,26 +302,21 @@ template void problem_checking_t::check_unscaled_solution( detail::problem_t& op_problem, rmm::device_uvector const& assignment) { - auto& d_variable_upper_bounds = op_problem.variable_upper_bounds; - auto& d_variable_lower_bounds = op_problem.variable_lower_bounds; - std::vector variable_upper_bounds(d_variable_upper_bounds.size()); - std::vector variable_lower_bounds(d_variable_lower_bounds.size()); + using f_t2 = typename type_2::type; + auto& d_variable_bounds = op_problem.variable_bounds; + std::vector variable_bounds(d_variable_bounds.size()); std::vector h_assignment(assignment.size()); - raft::copy(variable_upper_bounds.data(), - d_variable_upper_bounds.data(), - d_variable_upper_bounds.size(), - op_problem.handle_ptr->get_stream()); - raft::copy(variable_lower_bounds.data(), - d_variable_lower_bounds.data(), - d_variable_lower_bounds.size(), + raft::copy(variable_bounds.data(), + d_variable_bounds.data(), + d_variable_bounds.size(), op_problem.handle_ptr->get_stream()); raft::copy( h_assignment.data(), assignment.data(), assignment.size(), op_problem.handle_ptr->get_stream()); const f_t int_tol = op_problem.tolerances.integrality_tolerance; - for (size_t i = 0; i < variable_upper_bounds.size(); ++i) { - cuopt_assert(h_assignment[i] <= variable_upper_bounds[i] + int_tol, "Excess upper bound"); - cuopt_assert(h_assignment[i] >= variable_lower_bounds[i] - int_tol, "Excess lower bound"); + for (size_t i = 0; i < variable_bounds.size(); ++i) { + cuopt_assert(h_assignment[i] >= get_lower(variable_bounds[i]) - int_tol, "Excess lower bound"); + cuopt_assert(h_assignment[i] <= get_upper(variable_bounds[i]) + int_tol, "Excess upper bound"); } } diff --git a/cpp/src/linear_programming/utils.cuh b/cpp/src/linear_programming/utils.cuh index d4df2815b1..de107a66fb 100644 --- a/cpp/src/linear_programming/utils.cuh +++ b/cpp/src/linear_programming/utils.cuh @@ -77,13 +77,17 @@ struct a_sub_scalar_times_b { const f_t* scalar_; }; -template +template struct primal_projection { primal_projection(const f_t* step_size) : step_size_(step_size) {} - __device__ __forceinline__ thrust::tuple operator()( - f_t primal, f_t obj_coeff, f_t AtY, f_t lower, f_t upper) + __device__ __forceinline__ thrust::tuple operator()(f_t primal, + f_t obj_coeff, + f_t AtY, + f_t2 bounds) { + f_t lower = get_lower(bounds); + f_t upper = get_upper(bounds); f_t gradient = obj_coeff - AtY; f_t next = primal - (*step_size_ * gradient); next = raft::max(raft::min(next, upper), lower); @@ -129,13 +133,21 @@ struct a_divides_sqrt_b_bounded { }; template -struct clamp { +struct constraint_clamp { __device__ f_t operator()(f_t value, f_t lower, f_t upper) { return raft::min(raft::max(value, lower), upper); } }; +template +struct clamp { + __device__ f_t operator()(f_t value, f_t2 bounds) + { + return raft::min(raft::max(value, get_lower(bounds)), get_upper(bounds)); + } +}; + template struct combine_finite_abs_bounds { __device__ __host__ f_t operator()(f_t lower, f_t upper) @@ -177,32 +189,47 @@ struct violation { } }; -template +template struct max_violation { max_violation() {} - __device__ f_t operator()(const thrust::tuple& t) const + __device__ f_t operator()(const thrust::tuple& t) const { - const f_t value = thrust::get<0>(t); - const f_t lower = thrust::get<1>(t); - const f_t upper = thrust::get<2>(t); - f_t local_max = f_t(0.0); + const f_t value = thrust::get<0>(t); + const f_t2 bounds = thrust::get<1>(t); + const f_t lower = get_lower(bounds); + const f_t upper = get_upper(bounds); + f_t local_max = f_t(0.0); if (isfinite(lower)) { local_max = raft::max(local_max, -value); } if (isfinite(upper)) { local_max = raft::max(local_max, value); } return local_max; } }; -template +template +struct divide_check_zero { + __device__ f_t2 operator()(f_t2 bounds, f_t value) + { + if (value == f_t{0}) { + return f_t2{0, 0}; + } else { + return f_t2{get_lower(bounds) / value, get_upper(bounds) / value}; + } + } +}; + +template struct bound_value_gradient { - __device__ f_t operator()(f_t value, f_t lower, f_t upper) + __device__ f_t operator()(f_t value, f_t2 bounds) { + f_t lower = get_lower(bounds); + f_t upper = get_upper(bounds); if (value > f_t(0) && value < f_t(0)) { return 0; } return value > f_t(0) ? lower : upper; } }; template -struct bound_value_reduced_cost_product { +struct constraint_bound_value_reduced_cost_product { __device__ f_t operator()(f_t value, f_t lower, f_t upper) { f_t bound_value = f_t(0); @@ -218,6 +245,25 @@ struct bound_value_reduced_cost_product { } }; +template +struct bound_value_reduced_cost_product { + __device__ f_t operator()(f_t value, f_t2 variable_bounds) + { + f_t lower = get_lower(variable_bounds); + f_t upper = get_upper(variable_bounds); + f_t bound_value = f_t(0); + if (value > f_t(0)) { + // A positive reduced cost is associated with a binding lower bound. + bound_value = lower; + } else if (value < f_t(0)) { + // A negative reduced cost is associated with a binding upper bound. + bound_value = upper; + } + f_t val = isfinite(bound_value) ? value * bound_value : f_t(0); + return val; + } +}; + template struct copy_gradient_if_should_be_reduced_cost { __device__ f_t operator()(f_t value, f_t bound, f_t gradient) diff --git a/cpp/src/mip/CMakeLists.txt b/cpp/src/mip/CMakeLists.txt index 8f859e2d0d..7165b3bb56 100644 --- a/cpp/src/mip/CMakeLists.txt +++ b/cpp/src/mip/CMakeLists.txt @@ -36,19 +36,15 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/local_search/local_search.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/bounds_repair.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/constraint_prop.cu - ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/lb_bounds_repair.cu - ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/lb_constraint_prop.cu + ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/simple_rounding.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/feasibility_pump/feasibility_pump.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/line_segment_search/line_segment_search.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_presolve.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_update_data.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conditional_bound_strengthening.cu - ${CMAKE_CURRENT_SOURCE_DIR}/presolve/lb_probing_cache.cu - ${CMAKE_CURRENT_SOURCE_DIR}/presolve/load_balanced_bounds_presolve.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/multi_probe.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/probing_cache.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/trivial_presolve.cu - ${CMAKE_CURRENT_SOURCE_DIR}/problem/load_balanced_problem.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump_kernels.cu ) diff --git a/cpp/src/mip/diversity/assignment_hash_map.cu b/cpp/src/mip/diversity/assignment_hash_map.cu index 91ef05bd1f..24d0051b37 100644 --- a/cpp/src/mip/diversity/assignment_hash_map.cu +++ b/cpp/src/mip/diversity/assignment_hash_map.cu @@ -91,6 +91,7 @@ template size_t assignment_hash_map_t::hash_solution(solution_t& solution) { const int TPB = 1024; + fill_integer_assignment(solution); thrust::fill( solution.handle_ptr->get_thrust_policy(), reduction_buffer.begin(), reduction_buffer.end(), 0); diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu index b406d56a39..726eb5b41a 100644 --- a/cpp/src/mip/diversity/diversity_manager.cu +++ b/cpp/src/mip/diversity/diversity_manager.cu @@ -402,13 +402,13 @@ solution_t diversity_manager_t::run_solver() 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 + ls.resize_vectors(*problem_ptr, problem_ptr->handle_ptr); - ls.lb_constraint_prop.temp_problem.setup(*problem_ptr); - 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)); @@ -428,9 +428,6 @@ solution_t diversity_manager_t::run_solver() 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 = problem_ptr->lp_state; @@ -482,6 +479,7 @@ solution_t diversity_manager_t::run_solver() // 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(); } if (!fp_only_run) { @@ -505,7 +503,9 @@ solution_t diversity_manager_t::run_solver() } if (timer.check_time_limit()) { return population.best_feasible(); } + main_loop(); + return population.best_feasible(); }; diff --git a/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh b/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh index f38cc5759e..52347fc8b7 100644 --- a/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh +++ b/cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh @@ -86,8 +86,9 @@ class bound_prop_recombiner_t : public recombiner_t { 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]; + auto bounds = guiding_view.problem.variable_bounds[idx]; + f_t lb = get_lower(bounds); + f_t ub = get_upper(bounds); if (integer_equal(lb, ub, int_tol)) { cuopt_assert(false, "The var values must be different in A and B!"); } else if (isfinite(lb)) { diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump.cu b/cpp/src/mip/feasibility_jump/feasibility_jump.cu index 8e864dcf20..7520d1431e 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump.cu +++ b/cpp/src/mip/feasibility_jump/feasibility_jump.cu @@ -291,8 +291,8 @@ void fj_t::device_init(const rmm::cuda_stream_view& stream) cuopt_assert(var_idx < pb.is_binary_variable.size(), ""); if (pb.is_binary_variable[var_idx]) { - cuopt_assert(pb.variable_lower_bounds[var_idx] == 0 && - pb.variable_upper_bounds[var_idx] == 1, + cuopt_assert(get_lower(pb.variable_bounds[var_idx]) == 0 && + get_upper(pb.variable_bounds[var_idx]) == 1, "invalid bounds for binary variable"); } }); @@ -392,9 +392,9 @@ void fj_t::climber_init(i_t climber_idx, const rmm::cuda_stream_view& incumbent_assignment[var_idx] = round(incumbent_assignment[var_idx]); } // clamp to bounds + auto bounds = pb.variable_bounds[var_idx]; incumbent_assignment[var_idx] = - max(pb.variable_lower_bounds[var_idx], - min(pb.variable_upper_bounds[var_idx], incumbent_assignment[var_idx])); + max(get_lower(bounds), min(get_upper(bounds), incumbent_assignment[var_idx])); }); thrust::for_each( @@ -1102,6 +1102,7 @@ i_t fj_t::solve(solution_t& solution) 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(); diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu index 213d31bc5a..0d96cd4b2a 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu +++ b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu @@ -438,16 +438,15 @@ DI std::pair::move_score_info_t> compute_best_mtm( auto best_score_info = fj_t::move_score_info_t::invalid(); // fixed variables - if (fj.pb.integer_equal(fj.pb.variable_lower_bounds[var_idx], - fj.pb.variable_upper_bounds[var_idx])) { - return std::make_pair(fj.pb.variable_lower_bounds[var_idx], - fj_t::move_score_info_t::invalid()); + auto bounds = fj.pb.variable_bounds[var_idx]; + if (fj.pb.integer_equal(get_lower(bounds), get_upper(bounds))) { + return std::make_pair(get_lower(bounds), fj_t::move_score_info_t::invalid()); } f_t old_val = fj.incumbent_assignment[var_idx]; f_t obj_coeff = fj.pb.objective_coefficients[var_idx]; - f_t v_lb = fj.pb.variable_lower_bounds[var_idx]; - f_t v_ub = fj.pb.variable_upper_bounds[var_idx]; + f_t v_lb = get_lower(bounds); + f_t v_ub = get_upper(bounds); raft::random::PCGenerator rng(fj.settings->seed + *fj.iterations, 0, 0); cuopt_assert(isfinite(v_lb) || isfinite(v_ub), "unexpected free variable"); @@ -549,9 +548,9 @@ DI void update_jump_value(typename fj_t::climber_data_t::view_t fj, i_ 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, - ""); + cuopt_assert(get_lower(fj.pb.variable_bounds[var_idx]) == 0 && + get_upper(fj.pb.variable_bounds[var_idx]) == 1, + ""); cuopt_assert( fj.pb.check_variable_within_bounds(var_idx, fj.incumbent_assignment[var_idx] + delta), "Var not within bounds!"); @@ -566,8 +565,9 @@ DI void update_jump_value(typename fj_t::climber_data_t::view_t fj, i_ } else { 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, ""); + cuopt_assert(get_lower(fj.pb.variable_bounds[var_idx]) == 0 && + get_upper(fj.pb.variable_bounds[var_idx]) == 1, + ""); cuopt_assert( fj.pb.check_variable_within_bounds(var_idx, fj.incumbent_assignment[var_idx] + delta), "Var not within bounds!"); @@ -795,7 +795,8 @@ __global__ void update_assignment_kernel(typename fj_t::climber_data_t } } - i_t var_range = fj.pb.variable_upper_bounds[var_idx] - fj.pb.variable_lower_bounds[var_idx]; + auto bounds = fj.pb.variable_bounds[var_idx]; + i_t var_range = get_upper(bounds) - get_lower(bounds); 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; @@ -807,8 +808,8 @@ __global__ void update_assignment_kernel(typename fj_t::climber_data_t "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], + get_lower(fj.pb.variable_bounds[var_idx]), + get_upper(fj.pb.variable_bounds[var_idx]), fj.incumbent_assignment[var_idx], fj.jump_move_delta[var_idx], fj.incumbent_assignment[var_idx] + fj.jump_move_delta[var_idx], @@ -891,8 +892,9 @@ 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_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 bounds = fj.pb.variable_bounds[var_idx]; + f_t th_lower_delta = get_lower(bounds) - fj.incumbent_assignment[var_idx]; + f_t th_upper_delta = get_upper(bounds) - 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]; @@ -992,8 +994,9 @@ template DI f_t get_breakthrough_move(typename fj_t::climber_data_t::view_t fj, i_t var_idx) { f_t obj_coeff = fj.pb.objective_coefficients[var_idx]; - f_t v_lb = fj.pb.variable_lower_bounds[var_idx]; - f_t v_ub = fj.pb.variable_upper_bounds[var_idx]; + auto bounds = fj.pb.variable_bounds[var_idx]; + f_t v_lb = get_lower(bounds); + f_t v_ub = get_upper(bounds); cuopt_assert(isfinite(v_lb) || isfinite(v_ub), "unexpected free variable"); cuopt_assert(v_lb <= v_ub, "invalid bounds"); cuopt_assert(fj.pb.check_variable_within_bounds(var_idx, fj.incumbent_assignment[var_idx]), @@ -1217,8 +1220,8 @@ __device__ void compute_mtm_moves(typename fj_t::climber_data_t::view_ bool exclude_from_search = false; // "fixed" variables are to be excluded (as they cannot take any other value) - exclude_from_search |= fj.pb.integer_equal(fj.pb.variable_lower_bounds[var_idx], - fj.pb.variable_upper_bounds[var_idx]); + auto bounds = fj.pb.variable_bounds[var_idx]; + exclude_from_search |= fj.pb.integer_equal(get_lower(bounds), get_upper(bounds)); if (exclude_from_search) { if (threadIdx.x == 0) { @@ -1272,7 +1275,8 @@ __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]; + auto bounds = fj.pb.variable_bounds[var_idx]; + i_t var_range = get_upper(bounds) - get_lower(bounds); 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 && @@ -1319,16 +1323,16 @@ __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]; + auto bounds = fj.pb.variable_bounds[selected_var]; + i_t var_range = get_upper(bounds) - get_lower(bounds); double delta_rel_err = fabs(fj.jump_move_delta[selected_var]) / var_range * 100; DEVICE_LOG_INFO( "=---- 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], + get_lower(bounds), + get_upper(bounds), fj.pb.variable_types[selected_var] == var_t::INTEGER ? 'I' : 'C', fj.incumbent_assignment[selected_var], fj.jump_move_delta[selected_var], @@ -1553,9 +1557,10 @@ __global__ void handle_local_minimum_kernel(typename fj_t::climber_dat // 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], + auto bounds = fj.pb.variable_bounds[selected]; + delta = round_nearest(fj.incumbent_assignment[selected], + get_lower(bounds), + get_upper(bounds), fj.pb.tolerances.integrality_tolerance, rng) - fj.incumbent_assignment[selected]; @@ -1564,10 +1569,11 @@ __global__ void handle_local_minimum_kernel(typename fj_t::climber_dat if (FIRST_THREAD) { fj.jump_move_delta[selected] = delta; *fj.selected_var = selected; + auto bounds = fj.pb.variable_bounds[*fj.selected_var]; 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], + get_lower(bounds), + get_upper(bounds), fj.jump_move_delta[*fj.selected_var], fj.incumbent_assignment[*fj.selected_var]); } diff --git a/cpp/src/mip/feasibility_jump/load_balancing.cuh b/cpp/src/mip/feasibility_jump/load_balancing.cuh index b539687400..67af9c06ae 100644 --- a/cpp/src/mip/feasibility_jump/load_balancing.cuh +++ b/cpp/src/mip/feasibility_jump/load_balancing.cuh @@ -322,8 +322,9 @@ __global__ void load_balancing_compute_scores_binary( if (threadIdx.x == 0) { cuopt_assert(fj.incumbent_assignment[var_idx] == 0 || 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, ""); + cuopt_assert(get_lower(fj.pb.variable_bounds[var_idx]) == 0 && + get_upper(fj.pb.variable_bounds[var_idx]) == 1, + ""); cuopt_assert( fj.pb.check_variable_within_bounds(var_idx, fj.incumbent_assignment[var_idx] + delta), "Var not within bounds!"); @@ -400,8 +401,9 @@ __global__ void load_balancing_mtm_compute_candidates( auto rcp_cstr_coeff = fj.cstr_coeff_reciprocal[csr_offset]; f_t c_lb = fj.constraint_lower_bounds_csr[csr_offset]; f_t c_ub = fj.constraint_upper_bounds_csr[csr_offset]; - f_t v_lb = fj.pb.variable_lower_bounds[var_idx]; - f_t v_ub = fj.pb.variable_upper_bounds[var_idx]; + auto v_bnd = fj.pb.variable_bounds[var_idx]; + f_t v_lb = get_lower(v_bnd); + f_t v_ub = get_upper(v_bnd); cuopt_assert(c_lb == fj.pb.constraint_lower_bounds[cstr_idx], ""); cuopt_assert(c_ub == fj.pb.constraint_upper_bounds[cstr_idx], ""); @@ -512,8 +514,9 @@ __launch_bounds__(TPB_loadbalance, 16) __global__ cuopt_assert(cstr_idx >= 0 && cstr_idx < fj.pb.n_constraints, ""); } - f_t v_lb = fj.pb.variable_lower_bounds[var_idx]; - f_t v_ub = fj.pb.variable_upper_bounds[var_idx]; + auto v_bnd = fj.pb.variable_bounds[var_idx]; + f_t v_lb = get_lower(v_bnd); + f_t v_ub = get_upper(v_bnd); // candidate counts is usually very small (<4) thanks to early duplicate deletion in the // previous kernel rarely limits the thoroughput nor leads to noticeable imbalance 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 8286d8148e..f8b7d34fc5 100644 --- a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu +++ b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu @@ -46,7 +46,6 @@ feasibility_pump_t::feasibility_pump_t( fj_t& fj_, // fj_tree_t& fj_tree_, constraint_prop_t& constraint_prop_, - lb_constraint_prop_t& lb_constraint_prop_, line_segment_search_t& line_segment_search_, rmm::device_uvector& lp_optimal_solution_) : context(context_), @@ -55,7 +54,6 @@ feasibility_pump_t::feasibility_pump_t( line_segment_search(line_segment_search_), cycle_queue(*context.problem_ptr), constraint_prop(constraint_prop_), - lb_constraint_prop(lb_constraint_prop_), last_rounding(context.problem_ptr->n_variables, context.problem_ptr->handle_ptr->get_stream()), last_projection(context.problem_ptr->n_variables, context.problem_ptr->handle_ptr->get_stream()), @@ -147,11 +145,9 @@ bool feasibility_pump_t::linear_project_onto_polytope(solution_tvariable_upper_bounds, - solution.handle_ptr->get_stream()); - auto h_variable_lower_bounds = cuopt::host_copy(solution.problem_ptr->variable_lower_bounds, - solution.handle_ptr->get_stream()); + auto h_assignment = solution.get_host_assignment(); + auto h_variable_bounds = + cuopt::host_copy(solution.problem_ptr->variable_bounds, solution.handle_ptr->get_stream()); 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; @@ -164,23 +160,24 @@ bool feasibility_pump_t::linear_project_onto_polytope(solution_tinteger_equal(h_assignment[i], h_variable_upper_bounds[i])) { - obj_offset += h_variable_upper_bounds[i]; + auto h_var_bounds = h_variable_bounds[i]; + if (solution.problem_ptr->integer_equal(h_assignment[i], get_upper(h_var_bounds))) { + obj_offset += get_upper(h_var_bounds); // set the objective weight to -1, u - x obj_coefficients[i] = -1; - } else if (solution.problem_ptr->integer_equal(h_assignment[i], h_variable_lower_bounds[i])) { - obj_offset -= h_variable_lower_bounds[i]; + } else if (solution.problem_ptr->integer_equal(h_assignment[i], get_lower(h_var_bounds))) { + obj_offset -= get_lower(h_var_bounds); // set the objective weight to +1, x - l obj_coefficients[i] = 1; } else { // objective weight is 1 const f_t obj_weight = 1.; // the distance should always be positive - i_t var_id = h_variables.add_variable( - 0, - (h_variable_upper_bounds[i] - h_variable_lower_bounds[i]) + int_tol, - obj_weight, - var_t::CONTINUOUS); + i_t var_id = + h_variables.add_variable(0, + (get_upper(h_var_bounds) - get_lower(h_var_bounds)) + int_tol, + obj_weight, + var_t::CONTINUOUS); obj_coefficients.push_back(obj_weight); f_t dist_val = abs(h_assignment[i] - h_last_projection[i]); // if it is out of bounds, because of the approximation issues,or init issues @@ -442,8 +439,7 @@ void feasibility_pump_t::relax_general_integers(solution_t& orig_variable_types.resize(solution.problem_ptr->n_variables, solution.handle_ptr->get_stream()); auto var_types = make_span(solution.problem_ptr->variable_types); - auto var_lb = make_span(solution.problem_ptr->variable_lower_bounds); - auto var_ub = make_span(solution.problem_ptr->variable_upper_bounds); + auto var_bnds = make_span(solution.problem_ptr->variable_bounds); auto copy_types = make_span(orig_variable_types); raft::copy(orig_variable_types.data(), @@ -454,11 +450,11 @@ void feasibility_pump_t::relax_general_integers(solution_t& solution.handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), thrust::make_counting_iterator(solution.problem_ptr->n_variables), - [var_types, var_lb, var_ub, copy_types, pb = solution.problem_ptr->view()] __device__( - auto v_idx) { + [var_types, var_bnds, copy_types, pb = solution.problem_ptr->view()] __device__(auto v_idx) { auto orig_v_type = var_types[v_idx]; - auto lb = var_lb[v_idx]; - auto ub = var_ub[v_idx]; + auto var_bounds = var_bnds[v_idx]; + auto lb = get_lower(var_bounds); + auto ub = get_upper(var_bounds); bool var_binary = (pb.integer_equal(lb, 0) && pb.integer_equal(ub, 1)); auto copy_type = (orig_v_type == var_t::INTEGER) && var_binary ? var_t::INTEGER : var_t::CONTINUOUS; 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 60a249f893..2013e80f51 100644 --- a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cuh +++ b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cuh @@ -20,7 +20,6 @@ #include #include #include -#include #include #include @@ -118,7 +117,6 @@ class feasibility_pump_t { fj_t& fj, // fj_tree_t& fj_tree_, constraint_prop_t& constraint_prop_, - lb_constraint_prop_t& lb_constraint_prop_, line_segment_search_t& line_segment_search_, rmm::device_uvector& lp_optimal_solution_); @@ -153,7 +151,6 @@ class feasibility_pump_t { line_segment_search_t& line_segment_search; cycle_queue_t cycle_queue; constraint_prop_t& constraint_prop; - lb_constraint_prop_t& lb_constraint_prop; fp_config_t config; rmm::device_uvector last_rounding; rmm::device_uvector last_projection; diff --git a/cpp/src/mip/local_search/local_search.cu b/cpp/src/mip/local_search/local_search.cu index ae8a416fa9..ad444568fd 100644 --- a/cpp/src/mip/local_search/local_search.cu +++ b/cpp/src/mip/local_search/local_search.cu @@ -40,13 +40,11 @@ local_search_t::local_search_t(mip_solver_context_t& context fj(context), // fj_tree(fj), constraint_prop(context), - lb_constraint_prop(context), line_segment_search(fj, constraint_prop), fp(context, fj, // fj_tree, constraint_prop, - lb_constraint_prop, line_segment_search, lp_optimal_solution_), rng(cuopt::seed_generator::get_seed()), @@ -390,8 +388,6 @@ bool local_search_t::run_fp(solution_t& solution, 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) { @@ -441,8 +437,6 @@ bool local_search_t::run_fp(solution_t& solution, 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(); } diff --git a/cpp/src/mip/local_search/local_search.cuh b/cpp/src/mip/local_search/local_search.cuh index a2664e890b..bb95a8dc55 100644 --- a/cpp/src/mip/local_search/local_search.cuh +++ b/cpp/src/mip/local_search/local_search.cuh @@ -89,7 +89,6 @@ class local_search_t { fj_t fj; // fj_tree_t fj_tree; constraint_prop_t constraint_prop; - lb_constraint_prop_t lb_constraint_prop; line_segment_search_t line_segment_search; feasibility_pump_t fp; std::mt19937 rng; diff --git a/cpp/src/mip/local_search/rounding/bounds_repair.cu b/cpp/src/mip/local_search/rounding/bounds_repair.cu index a63c6c8669..e77702c38f 100644 --- a/cpp/src/mip/local_search/rounding/bounds_repair.cu +++ b/cpp/src/mip/local_search/rounding/bounds_repair.cu @@ -152,10 +152,12 @@ i_t bounds_repair_t::compute_best_shift(problem_t& problem, shift_amount = (down_vio / var_coeff); } if (shift_amount != 0.) { - f_t var_lb = pb_v.variable_lower_bounds[var_idx]; - f_t var_ub = pb_v.variable_upper_bounds[var_idx]; - f_t o_var_lb = o_pb_v.variable_lower_bounds[var_idx]; - f_t o_var_ub = o_pb_v.variable_upper_bounds[var_idx]; + auto var_bnd = pb_v.variable_bounds[var_idx]; + auto o_var_bnd = o_pb_v.variable_bounds[var_idx]; + f_t var_lb = get_lower(var_bnd); + f_t var_ub = get_upper(var_bnd); + f_t o_var_lb = get_lower(o_var_bnd); + f_t o_var_ub = get_upper(o_var_bnd); cuopt_assert(var_lb + pb_v.tolerances.integrality_tolerance >= o_var_lb, ""); cuopt_assert(o_var_ub + pb_v.tolerances.integrality_tolerance >= var_ub, ""); // round the shift amount of integer @@ -211,8 +213,9 @@ __global__ void compute_damages_kernel(typename problem_t::view_t prob { i_t var_idx = candidates.variable_index[blockIdx.x]; f_t shift_amount = candidates.bound_shift[blockIdx.x]; - f_t v_lb = problem.variable_lower_bounds[var_idx]; - f_t v_ub = problem.variable_upper_bounds[var_idx]; + auto v_bnd = problem.variable_bounds[var_idx]; + f_t v_lb = get_lower(v_bnd); + f_t v_ub = get_upper(v_bnd); f_t th_damage = 0.; i_t n_infeasible_cstr_delta = 0; auto [offset_begin, offset_end] = problem.reverse_range_for_var(var_idx); @@ -348,37 +351,37 @@ void bounds_repair_t::apply_move(problem_t& problem, problem_t& original_problem, i_t move_idx) { - run_device_lambda( - handle_ptr->get_stream(), - [move_idx, - candidates = candidates.view(), - problem = problem.view(), - original_problem = original_problem.view()] __device__() { - i_t var_idx = candidates.variable_index[move_idx]; - f_t shift_value = candidates.bound_shift[move_idx]; - DEVICE_LOG_TRACE("Applying move on var %d with shift %f lb %f ub %f o_lb %f o_ub %f \n", - var_idx, - shift_value, - problem.variable_lower_bounds[var_idx], - problem.variable_upper_bounds[var_idx], - original_problem.variable_lower_bounds[var_idx], - original_problem.variable_upper_bounds[var_idx]); - if (problem.integer_equal(problem.variable_lower_bounds[var_idx], - problem.variable_upper_bounds[var_idx])) { - *candidates.at_least_one_singleton_moved = 1; - } + run_device_lambda(handle_ptr->get_stream(), + [move_idx, + candidates = candidates.view(), + problem = problem.view(), + original_problem = original_problem.view()] __device__() { + i_t var_idx = candidates.variable_index[move_idx]; + f_t shift_value = candidates.bound_shift[move_idx]; + auto bounds = problem.variable_bounds[var_idx]; + DEVICE_LOG_TRACE( + "Applying move on var %d with shift %f lb %f ub %f o_lb %f o_ub %f \n", + var_idx, + shift_value, + get_lower(bounds), + get_upper(bounds), + get_lower(original_problem.variable_bounds[var_idx]), + get_upper(original_problem.variable_bounds[var_idx])); + if (problem.integer_equal(get_lower(bounds), get_upper(bounds))) { + *candidates.at_least_one_singleton_moved = 1; + } - problem.variable_lower_bounds[var_idx] += shift_value; - problem.variable_upper_bounds[var_idx] += shift_value; - cuopt_assert( - original_problem.variable_lower_bounds[var_idx] <= - problem.variable_lower_bounds[var_idx] + problem.tolerances.integrality_tolerance, - ""); - cuopt_assert(original_problem.variable_upper_bounds[var_idx] + - problem.tolerances.integrality_tolerance >= - problem.variable_upper_bounds[var_idx], - ""); - }); + get_lower(bounds) += shift_value; + get_upper(bounds) += shift_value; + problem.variable_bounds[var_idx] = bounds; + cuopt_assert(get_lower(original_problem.variable_bounds[var_idx]) <= + get_lower(bounds) + problem.tolerances.integrality_tolerance, + ""); + cuopt_assert(get_upper(original_problem.variable_bounds[var_idx]) + + problem.tolerances.integrality_tolerance >= + get_upper(bounds), + ""); + }); } template diff --git a/cpp/src/mip/local_search/rounding/bounds_repair.cuh b/cpp/src/mip/local_search/rounding/bounds_repair.cuh index 2a57d06600..8cc392c0b9 100644 --- a/cpp/src/mip/local_search/rounding/bounds_repair.cuh +++ b/cpp/src/mip/local_search/rounding/bounds_repair.cuh @@ -43,17 +43,25 @@ struct bounds_t { } void update_from(const problem_t& pb, const raft::handle_t* handle_ptr) { - cuopt_assert(lb.size() == pb.variable_lower_bounds.size(), ""); - cuopt_assert(ub.size() == pb.variable_upper_bounds.size(), ""); - raft::copy(lb.data(), pb.variable_lower_bounds.data(), lb.size(), handle_ptr->get_stream()); - raft::copy(ub.data(), pb.variable_upper_bounds.data(), ub.size(), handle_ptr->get_stream()); + cuopt_assert(lb.size() == pb.variable_bounds.size(), ""); + cuopt_assert(ub.size() == pb.variable_bounds.size(), ""); + thrust::transform( + handle_ptr->get_thrust_policy(), + pb.variable_bounds.begin(), + pb.variable_bounds.end(), + thrust::make_zip_iterator(thrust::make_tuple(lb.begin(), ub.begin())), + [] __device__(auto i) { return thrust::make_tuple(get_lower(i), get_upper(i)); }); }; void update_to(problem_t& pb, const raft::handle_t* handle_ptr) { - cuopt_assert(lb.size() == pb.variable_lower_bounds.size(), ""); - cuopt_assert(ub.size() == pb.variable_upper_bounds.size(), ""); - raft::copy(pb.variable_lower_bounds.data(), lb.data(), lb.size(), handle_ptr->get_stream()); - raft::copy(pb.variable_upper_bounds.data(), ub.data(), ub.size(), handle_ptr->get_stream()); + cuopt_assert(lb.size() == pb.variable_bounds.size(), ""); + cuopt_assert(ub.size() == pb.variable_bounds.size(), ""); + using f_t2 = typename type_2::type; + thrust::transform(handle_ptr->get_thrust_policy(), + thrust::make_zip_iterator(thrust::make_tuple(lb.begin(), ub.begin())), + thrust::make_zip_iterator(thrust::make_tuple(lb.end(), ub.end())), + pb.variable_bounds.begin(), + [] __device__(auto i) { return f_t2{thrust::get<0>(i), thrust::get<1>(i)}; }); }; rmm::device_uvector lb; rmm::device_uvector ub; diff --git a/cpp/src/mip/local_search/rounding/constraint_prop.cu b/cpp/src/mip/local_search/rounding/constraint_prop.cu index 6e6a5deb37..61e8e08675 100644 --- a/cpp/src/mip/local_search/rounding/constraint_prop.cu +++ b/cpp/src/mip/local_search/rounding/constraint_prop.cu @@ -130,8 +130,9 @@ __global__ void compute_implied_slack_consumption_per_var( i_t var_offset = pb.reverse_offsets[var_idx]; i_t var_degree = pb.reverse_offsets[var_idx + 1] - var_offset; f_t th_var_implied_slack_consumption = 0.; - f_t lb = pb.variable_lower_bounds[var_idx]; - f_t ub = pb.variable_upper_bounds[var_idx]; + auto var_bnd = pb.variable_bounds[var_idx]; + f_t lb = get_lower(var_bnd); + f_t ub = get_upper(var_bnd); for (i_t i = threadIdx.x; i < var_degree; i += blockDim.x) { auto a = pb.reverse_coefficients[var_offset + i]; auto cnst_idx = pb.reverse_constraints[var_offset + i]; @@ -206,25 +207,26 @@ void constraint_prop_t::sort_by_interval_and_frac(solution_t // we can't call this function when the problem is ii. it causes false offset computations // TODO add assert that the problem is not ii auto assgn = make_span(sol.assignment); - thrust::stable_sort(sol.handle_ptr->get_thrust_policy(), - vars.begin(), - vars.end(), - [lb = sol.problem_ptr->variable_lower_bounds.data(), - ub = sol.problem_ptr->variable_upper_bounds.data(), - assgn] __device__(i_t v_idx_1, i_t v_idx_2) { - f_t bounds_interval_1 = ub[v_idx_1] - lb[v_idx_1]; - f_t bounds_interval_2 = ub[v_idx_2] - lb[v_idx_2]; - // if bounds interval are equal (binary and ternary) check fraction - // if both bounds intervals are greater than 2. then do fraction - if ((bounds_interval_1 == bounds_interval_2) || - (bounds_interval_1 > 2 && bounds_interval_2 > 2)) { - f_t frac_1 = get_fractionality_of_val(assgn[v_idx_1]); - f_t frac_2 = get_fractionality_of_val(assgn[v_idx_2]); - return frac_1 < frac_2; - } else { - return bounds_interval_1 < bounds_interval_2; - } - }); + thrust::stable_sort( + sol.handle_ptr->get_thrust_policy(), + vars.begin(), + vars.end(), + [bnds = sol.problem_ptr->variable_bounds.data(), assgn] __device__(i_t v_idx_1, i_t v_idx_2) { + auto bnd_1 = bnds[v_idx_1]; + auto bnd_2 = bnds[v_idx_2]; + f_t bounds_interval_1 = get_upper(bnd_1) - get_lower(bnd_1); + f_t bounds_interval_2 = get_upper(bnd_2) - get_lower(bnd_2); + // if bounds interval are equal (binary and ternary) check fraction + // if both bounds intervals are greater than 2. then do fraction + if ((bounds_interval_1 == bounds_interval_2) || + (bounds_interval_1 > 2 && bounds_interval_2 > 2)) { + f_t frac_1 = get_fractionality_of_val(assgn[v_idx_1]); + f_t frac_2 = get_fractionality_of_val(assgn[v_idx_2]); + return frac_1 < frac_2; + } else { + return bounds_interval_1 < bounds_interval_2; + } + }); // now do the suffling, for that we need to assign some random values to rnd array // we will sort this rnd array and the vars in subsections, so that each subsection will be // shuffled in total we will have 3(binary, ternary and rest) x 7 intervals = 21 subsections. @@ -237,15 +239,16 @@ void constraint_prop_t::sort_by_interval_and_frac(solution_t thrust::for_each(sol.handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), thrust::make_counting_iterator((i_t)vars.size() - 1), - [lb = make_span(sol.problem_ptr->variable_lower_bounds), - ub = make_span(sol.problem_ptr->variable_upper_bounds), + [bnds = make_span(sol.problem_ptr->variable_bounds), offsets = make_span(subsection_offsets), vars, assgn] __device__(i_t idx) { i_t var_1 = vars[idx]; i_t var_2 = vars[idx + 1]; - f_t bounds_interval_1 = ub[var_1] - lb[var_1]; - f_t bounds_interval_2 = ub[var_2] - lb[var_2]; + auto bnd_1 = bnds[var_1]; + auto bnd_2 = bnds[var_2]; + f_t bounds_interval_1 = get_upper(bnd_1) - get_lower(bnd_1); + f_t bounds_interval_2 = get_upper(bnd_2) - get_lower(bnd_2); f_t frac_1 = get_fractionality_of_val(assgn[var_1]); f_t frac_2 = get_fractionality_of_val(assgn[var_2]); if (bounds_interval_1 == 1 && bounds_interval_2 == 1) { @@ -390,24 +393,22 @@ void constraint_prop_t::collapse_crossing_bounds(problem_t& problem_t& orig_problem, const raft::handle_t* handle_ptr) { - auto lb = make_span(problem.variable_lower_bounds); - auto ub = make_span(problem.variable_upper_bounds); - auto original_lb = make_span(orig_problem.variable_lower_bounds); - auto original_ub = make_span(orig_problem.variable_upper_bounds); + auto v_bnds = make_span(problem.variable_bounds); + auto original_v_bnds = make_span(orig_problem.variable_bounds); thrust::for_each( handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), - thrust::make_counting_iterator((i_t)lb.size()), - [lb, - ub, - original_lb, - original_ub, + thrust::make_counting_iterator((i_t)v_bnds.size()), + [v_bnds, + original_v_bnds, variable_types = make_span(problem.variable_types), int_tol = problem.tolerances.integrality_tolerance] __device__(i_t idx) { - auto v_lb = lb[idx]; - auto v_ub = ub[idx]; - auto o_lb = original_lb[idx]; - auto o_ub = original_ub[idx]; + auto v_bnd = v_bnds[idx]; + auto ov_bnd = original_v_bnds[idx]; + auto v_lb = get_lower(v_bnd); + auto v_ub = get_upper(v_bnd); + auto o_lb = get_lower(ov_bnd); + auto o_ub = get_upper(ov_bnd); if (v_lb > v_ub) { f_t val_to_collapse; if (variable_types[idx] == var_t::INTEGER) { @@ -422,8 +423,7 @@ void constraint_prop_t::collapse_crossing_bounds(problem_t& cuopt_assert(o_lb - int_tol <= val_to_collapse && val_to_collapse <= o_ub + int_tol, "Out of original bounds!"); - lb[idx] = val_to_collapse; - ub[idx] = val_to_collapse; + v_bnds[idx] = typename type_2::type{val_to_collapse, val_to_collapse}; } }); } @@ -431,51 +431,42 @@ void constraint_prop_t::collapse_crossing_bounds(problem_t& template void constraint_prop_t::set_bounds_on_fixed_vars(solution_t& sol) { - auto assgn = make_span(sol.assignment); - auto lb = make_span(sol.problem_ptr->variable_lower_bounds); - auto ub = make_span(sol.problem_ptr->variable_upper_bounds); + auto assgn = make_span(sol.assignment); + auto var_bounds = make_span(sol.problem_ptr->variable_bounds); thrust::for_each(sol.handle_ptr->get_thrust_policy(), sol.problem_ptr->integer_indices.begin(), sol.problem_ptr->integer_indices.end(), - [pb = sol.problem_ptr->view(), assgn, lb, ub] __device__(i_t idx) { + [pb = sol.problem_ptr->view(), assgn, var_bounds] __device__(i_t idx) { auto var_val = assgn[idx]; if (pb.is_integer(var_val)) { - lb[idx] = var_val; - ub[idx] = var_val; + var_bounds[idx] = typename type_2::type{var_val, var_val}; } }); } -template +template struct is_bound_fixed_t { // This functor should be called only on integer variables f_t eps; - raft::device_span lb; - raft::device_span ub; - raft::device_span original_lb; - raft::device_span original_ub; + raft::device_span::type> bnd; + raft::device_span::type> original_bnd; raft::device_span assignment; is_bound_fixed_t(f_t eps_, - raft::device_span lb_, - raft::device_span ub_, - raft::device_span original_lb_, - raft::device_span original_ub_, + raft::device_span bnd_, + raft::device_span original_bnd_, raft::device_span assignment_) - : eps(eps_), - lb(lb_), - ub(ub_), - original_lb(original_lb_), - original_ub(original_ub_), - assignment(assignment_) + : eps(eps_), bnd(bnd_), original_bnd(original_bnd_), assignment(assignment_) { } HDI bool operator()(i_t idx) { - auto v_lb = lb[idx]; - auto v_ub = ub[idx]; - auto o_lb = original_lb[idx]; - auto o_ub = original_ub[idx]; + auto v_bnd = bnd[idx]; + auto v_lb = get_lower(v_bnd); + auto v_ub = get_upper(v_bnd); + auto ov_bnd = original_bnd[idx]; + auto o_lb = get_lower(ov_bnd); + auto o_ub = get_upper(ov_bnd); bool is_singleton = round_val_on_singleton_and_crossing(assignment[idx], v_lb, v_ub, o_lb, o_ub); return is_singleton; @@ -518,6 +509,40 @@ struct greater_than_threshold_t { __host__ __device__ bool operator()(const i_t& x) const { return assignment[x] > threshold; } }; +template +void constraint_prop_t::copy_bounds( + rmm::device_uvector::type>& output_bounds, + const rmm::device_uvector& input_lb, + const rmm::device_uvector& input_ub, + const raft::handle_t* handle_ptr) +{ + thrust::transform( + handle_ptr->get_thrust_policy(), + thrust::make_zip_iterator(thrust::make_tuple(input_lb.begin(), input_ub.begin())), + thrust::make_zip_iterator(thrust::make_tuple(input_lb.end(), input_ub.end())), + output_bounds.begin(), + [] __device__(auto bounds) { + return typename type_2::type{thrust::get<0>(bounds), thrust::get<1>(bounds)}; + }); +} + +template +void constraint_prop_t::copy_bounds( + rmm::device_uvector& output_lb, + rmm::device_uvector& output_ub, + const rmm::device_uvector::type>& input_bounds, + const raft::handle_t* handle_ptr) +{ + thrust::transform( + handle_ptr->get_thrust_policy(), + input_bounds.begin(), + input_bounds.end(), + thrust::make_zip_iterator(thrust::make_tuple(output_lb.begin(), output_ub.begin())), + [] __device__(auto bounds) { + return thrust::make_tuple(get_lower(bounds), get_upper(bounds)); + }); +} + template void constraint_prop_t::copy_bounds(rmm::device_uvector& output_lb, rmm::device_uvector& output_ub, @@ -548,38 +573,35 @@ void constraint_prop_t::copy_bounds(rmm::device_uvector& output_l template void constraint_prop_t::save_bounds(solution_t& sol) { - copy_bounds(lb_restore, - ub_restore, - assignment_restore, - sol.problem_ptr->variable_lower_bounds, - sol.problem_ptr->variable_upper_bounds, - sol.assignment, - sol.handle_ptr); + copy_bounds(lb_restore, ub_restore, sol.problem_ptr->variable_bounds, sol.handle_ptr); + raft::copy(assignment_restore.data(), + sol.assignment.data(), + sol.assignment.size(), + sol.handle_ptr->get_stream()); } template void constraint_prop_t::restore_bounds(solution_t& sol) { - copy_bounds(sol.problem_ptr->variable_lower_bounds, - sol.problem_ptr->variable_upper_bounds, - sol.assignment, - lb_restore, - ub_restore, - assignment_restore, - sol.handle_ptr); + copy_bounds(sol.problem_ptr->variable_bounds, lb_restore, ub_restore, sol.handle_ptr); + raft::copy(sol.assignment.data(), + assignment_restore.data(), + assignment_restore.size(), + sol.handle_ptr->get_stream()); } template void constraint_prop_t::restore_original_bounds(solution_t& sol, solution_t& orig_sol) { - copy_bounds(sol.problem_ptr->variable_lower_bounds, - sol.problem_ptr->variable_upper_bounds, - sol.assignment, - orig_sol.problem_ptr->variable_lower_bounds, - orig_sol.problem_ptr->variable_upper_bounds, - orig_sol.assignment, - orig_sol.handle_ptr); + raft::copy(sol.problem_ptr->variable_bounds.data(), + orig_sol.problem_ptr->variable_bounds.data(), + orig_sol.problem_ptr->variable_bounds.size(), + orig_sol.handle_ptr->get_stream()); + raft::copy(sol.assignment.data(), + orig_sol.assignment.data(), + orig_sol.assignment.size(), + orig_sol.handle_ptr->get_stream()); } template @@ -635,6 +657,18 @@ thrust::pair constraint_prop_t::generate_double_probing_pair return thrust::make_pair(first_probe, second_probe); } +template +bool test_var_out_of_bounds(const solution_t& orig_sol, + i_t unset_var_idx, + f_t probe, + f_t int_tol, + const raft::handle_t* handle_ptr) +{ + auto var_bnd = + orig_sol.problem_ptr->variable_bounds.element(unset_var_idx, handle_ptr->get_stream()); + return (get_lower(var_bnd) <= probe + int_tol) && (probe - int_tol <= get_upper(var_bnd)); +} + template std::tuple, std::vector, std::vector> constraint_prop_t::generate_bulk_rounding_vector( @@ -660,16 +694,12 @@ constraint_prop_t::generate_bulk_rounding_vector( 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->variable_lower_bounds.element( - unset_var_idx, sol.handle_ptr->get_stream()) <= first_probe + int_tol && - first_probe - int_tol <= orig_sol.problem_ptr->variable_upper_bounds.element( - unset_var_idx, sol.handle_ptr->get_stream()), - "Variable out of original bounds!"); - 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!"); + cuopt_assert( + test_var_out_of_bounds(orig_sol, unset_var_idx, first_probe, int_tol, sol.handle_ptr), + "Variable out of original bounds!"); + cuopt_assert( + test_var_out_of_bounds(orig_sol, unset_var_idx, second_probe, int_tol, sol.handle_ptr), + "Variable out of original bounds!"); 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), "Probing value must be an integer"); @@ -687,16 +717,12 @@ constraint_prop_t::generate_bulk_rounding_vector( 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!"); - 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!"); + cuopt_assert( + test_var_out_of_bounds(orig_sol, unset_var_idx, val_to_round, int_tol, sol.handle_ptr), + "Variable out of original bounds!"); + cuopt_assert( + test_var_out_of_bounds(orig_sol, unset_var_idx, second_probe, int_tol, sol.handle_ptr), + "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; @@ -718,18 +744,8 @@ 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() == multi_probe.host_lb.size(), - "size of variable lower bound mismatch"); - 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() == multi_probe.host_ub.size(), - "size of variable upper bound mismatch"); - 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()); + std::tie(multi_probe.host_lb, multi_probe.host_ub) = + extract_host_bounds(sol.problem_ptr->variable_bounds, sol.handle_ptr); } template @@ -742,11 +758,10 @@ void constraint_prop_t::restore_original_bounds_on_unfixed( thrust::make_counting_iterator(0), thrust::make_counting_iterator(problem.n_variables), [p_v = problem.view(), op_v = original_problem.view()] __device__(i_t var_idx) { - if (!p_v.integer_equal(p_v.variable_lower_bounds[var_idx], - p_v.variable_upper_bounds[var_idx]) || + auto p_v_var_bnd = p_v.variable_bounds[var_idx]; + if (!p_v.integer_equal(get_lower(p_v_var_bnd), get_upper(p_v_var_bnd)) || !p_v.is_integer_var(var_idx)) { - p_v.variable_lower_bounds[var_idx] = op_v.variable_lower_bounds[var_idx]; - p_v.variable_upper_bounds[var_idx] = op_v.variable_upper_bounds[var_idx]; + p_v.variable_bounds[var_idx] = op_v.variable_bounds[var_idx]; } }); } @@ -947,15 +962,14 @@ bool constraint_prop_t::find_integer( rounding_ii = false; n_iter_in_recovery = 0; // during repair procedure some variables might be collapsed - auto iter = thrust::stable_partition( - sol.handle_ptr->get_thrust_policy(), - unset_vars.begin() + set_count, - unset_vars.end(), - is_bound_fixed_t{orig_sol.problem_ptr->tolerances.integrality_tolerance, - make_span(sol.problem_ptr->variable_lower_bounds), - make_span(sol.problem_ptr->variable_upper_bounds), - make_span(orig_sol.problem_ptr->variable_lower_bounds), - make_span(orig_sol.problem_ptr->variable_upper_bounds), + auto iter = + thrust::stable_partition(sol.handle_ptr->get_thrust_policy(), + unset_vars.begin() + set_count, + unset_vars.end(), + is_bound_fixed_t::type>{ + orig_sol.problem_ptr->tolerances.integrality_tolerance, + make_span(sol.problem_ptr->variable_bounds), + make_span(orig_sol.problem_ptr->variable_bounds), make_span(sol.assignment)}); i_t n_fixed_vars = (iter - (unset_vars.begin() + set_count)); CUOPT_LOG_TRACE("After repair procedure, number of additional fixed vars %d", n_fixed_vars); @@ -983,9 +997,7 @@ 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 - multi_probe.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_bounds)); } CUOPT_LOG_DEBUG( "Bounds propagation rounding end: ii constraint count first buffer %d, second buffer %d", @@ -1096,10 +1108,10 @@ std::tuple constraint_prop_t::probing_values( "probing value out of bounds"); return std::make_tuple(first_round_val, var_val, second_round_val); } else { - auto orig_v_lb = - orig_sol.problem_ptr->variable_lower_bounds.element(idx, sol.handle_ptr->get_stream()); - auto orig_v_ub = - orig_sol.problem_ptr->variable_upper_bounds.element(idx, sol.handle_ptr->get_stream()); + auto orig_v_bnd = + orig_sol.problem_ptr->variable_bounds.element(idx, sol.handle_ptr->get_stream()); + auto orig_v_lb = get_lower(orig_v_bnd); + auto orig_v_ub = get_upper(orig_v_bnd); cuopt_assert(v_lb >= orig_v_lb, "Current lb should be greater than original lb"); cuopt_assert(v_ub <= orig_v_ub, "Current ub should be smaller than original ub"); v_lb = std::max(v_lb, orig_v_lb); @@ -1137,16 +1149,14 @@ bool constraint_prop_t::handle_fixed_vars( 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(), - unset_vars.begin() + set_count, - unset_vars.end(), - is_bound_fixed_t{int_tol, - make_span(sol.problem_ptr->variable_lower_bounds), - make_span(sol.problem_ptr->variable_upper_bounds), - make_span(original_problem->variable_lower_bounds), - make_span(original_problem->variable_upper_bounds), - make_span(sol.assignment)}); + auto iter = thrust::stable_partition(sol.handle_ptr->get_thrust_policy(), + unset_vars.begin() + set_count, + unset_vars.end(), + is_bound_fixed_t::type>{ + int_tol, + make_span(sol.problem_ptr->variable_bounds), + make_span(original_problem->variable_bounds), + make_span(sol.assignment)}); i_t n_fixed_vars = (iter - (unset_vars.begin() + set_count)); cuopt_assert(n_fixed_vars >= std::get<0>(var_probe_vals).size(), "Error in number of vars fixed!"); diff --git a/cpp/src/mip/local_search/rounding/constraint_prop.cuh b/cpp/src/mip/local_search/rounding/constraint_prop.cuh index cb9cbbd00b..3b01da2749 100644 --- a/cpp/src/mip/local_search/rounding/constraint_prop.cuh +++ b/cpp/src/mip/local_search/rounding/constraint_prop.cuh @@ -95,6 +95,17 @@ 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::type>& input_bounds, + const raft::handle_t* handle_ptr); + + void copy_bounds(rmm::device_uvector::type>& output_bounds, + 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, const rmm::device_uvector& input_lb, diff --git a/cpp/src/mip/local_search/rounding/simple_rounding_kernels.cuh b/cpp/src/mip/local_search/rounding/simple_rounding_kernels.cuh index 9c1ac7f3b1..1c0103110c 100644 --- a/cpp/src/mip/local_search/rounding/simple_rounding_kernels.cuh +++ b/cpp/src/mip/local_search/rounding/simple_rounding_kernels.cuh @@ -36,8 +36,9 @@ __global__ void nearest_rounding_kernel(typename solution_t::view_t so f_t curr_val = solution.assignment[var_id]; if (solution.problem.is_integer(curr_val)) { return; } const f_t int_tol = solution.problem.tolerances.integrality_tolerance; - f_t lb = solution.problem.variable_lower_bounds[var_id]; - f_t ub = solution.problem.variable_upper_bounds[var_id]; + auto var_bnd = solution.problem.variable_bounds[var_id]; + f_t lb = get_lower(var_bnd); + f_t ub = get_upper(var_bnd); f_t nearest_val = round_nearest(curr_val, lb, ub, int_tol, rng); solution.assignment[var_id] = nearest_val; } diff --git a/cpp/src/mip/presolve/bounds_presolve.cu b/cpp/src/mip/presolve/bounds_presolve.cu index 72440cd9a8..35d2ee0821 100644 --- a/cpp/src/mip/presolve/bounds_presolve.cu +++ b/cpp/src/mip/presolve/bounds_presolve.cu @@ -212,14 +212,7 @@ void bound_presolve_t::calculate_activity_on_problem_bounds(problem_t< { auto& handle_ptr = pb.handle_ptr; upd.init_changed_constraints(handle_ptr); - cuopt_assert(upd.lb.size() == pb.variable_lower_bounds.size(), - "size of variable lower bound mismatch"); - raft::copy( - upd.lb.data(), pb.variable_lower_bounds.data(), upd.lb.size(), handle_ptr->get_stream()); - cuopt_assert(upd.ub.size() == pb.variable_upper_bounds.size(), - "size of variable upper bound mismatch"); - raft::copy( - upd.ub.data(), pb.variable_upper_bounds.data(), upd.ub.size(), handle_ptr->get_stream()); + copy_input_bounds(pb); calculate_activity(pb); } @@ -228,14 +221,15 @@ void bound_presolve_t::copy_input_bounds(problem_t& pb) { auto& handle_ptr = pb.handle_ptr; - cuopt_assert(upd.lb.size() == pb.variable_lower_bounds.size(), - "size of variable lower bound mismatch"); - raft::copy( - upd.lb.data(), pb.variable_lower_bounds.data(), upd.lb.size(), handle_ptr->get_stream()); - cuopt_assert(upd.ub.size() == pb.variable_upper_bounds.size(), - "size of variable upper bound mismatch"); - raft::copy( - upd.ub.data(), pb.variable_upper_bounds.data(), upd.ub.size(), handle_ptr->get_stream()); + cuopt_assert(upd.lb.size() == pb.variable_bounds.size(), "size of variable lower bound mismatch"); + cuopt_assert(upd.ub.size() == pb.variable_bounds.size(), "size of variable upper bound mismatch"); + + thrust::transform( + handle_ptr->get_thrust_policy(), + pb.variable_bounds.begin(), + pb.variable_bounds.end(), + thrust::make_zip_iterator(thrust::make_tuple(upd.lb.begin(), upd.ub.begin())), + [] __device__(auto i) { return thrust::make_tuple(get_lower(i), get_upper(i)); }); } template @@ -271,30 +265,11 @@ termination_criterion_t bound_presolve_t::solve( } template -termination_criterion_t bound_presolve_t::solve(problem_t& pb, - raft::device_span input_lb, - raft::device_span input_ub) +termination_criterion_t bound_presolve_t::solve(problem_t& pb) { timer_t timer(settings.time_limit); auto& handle_ptr = pb.handle_ptr; - if (input_lb.size() == 0) { - cuopt_assert(upd.lb.size() == pb.variable_lower_bounds.size(), - "size of variable lower bound mismatch"); - raft::copy( - upd.lb.data(), pb.variable_lower_bounds.data(), upd.lb.size(), handle_ptr->get_stream()); - } else { - cuopt_assert(input_lb.size() == upd.lb.size(), "size of variable lower bound mismatch"); - raft::copy(upd.lb.data(), input_lb.data(), input_lb.size(), handle_ptr->get_stream()); - } - if (input_ub.size() == 0) { - cuopt_assert(upd.ub.size() == pb.variable_upper_bounds.size(), - "size of variable upper bound mismatch"); - raft::copy( - upd.ub.data(), pb.variable_upper_bounds.data(), upd.ub.size(), handle_ptr->get_stream()); - } else { - cuopt_assert(input_ub.size() == upd.ub.size(), "size of variable lower bound mismatch"); - raft::copy(upd.ub.data(), input_ub.data(), upd.ub.size(), handle_ptr->get_stream()); - } + copy_input_bounds(pb); return bound_update_loop(pb, timer); } @@ -329,9 +304,7 @@ bool bound_presolve_t::calculate_infeasible_redundant_constraints(prob template void bound_presolve_t::set_updated_bounds(problem_t& pb) { - set_updated_bounds(pb.handle_ptr, - cuopt::make_span(pb.variable_lower_bounds), - cuopt::make_span(pb.variable_upper_bounds)); + set_updated_bounds(pb.handle_ptr, cuopt::make_span(pb.variable_bounds)); pb.compute_n_integer_vars(); pb.compute_binary_var_table(); } @@ -347,6 +320,21 @@ void bound_presolve_t::set_updated_bounds(const raft::handle_t* handle raft::copy(output_ub.data(), upd.ub.data(), upd.ub.size(), handle_ptr->get_stream()); } +template +void bound_presolve_t::set_updated_bounds( + const raft::handle_t* handle_ptr, raft::device_span::type> output_bounds) +{ + cuopt_assert(upd.ub.size() == output_bounds.size(), "size of variable upper bound mismatch"); + cuopt_assert(upd.lb.size() == output_bounds.size(), "size of variable lower bound mismatch"); + thrust::transform(handle_ptr->get_thrust_policy(), + thrust::make_zip_iterator(thrust::make_tuple(upd.lb.begin(), upd.ub.begin())), + thrust::make_zip_iterator(thrust::make_tuple(upd.lb.end(), upd.ub.end())), + output_bounds.begin(), + [] __device__(auto i) { + return typename type_2::type{thrust::get<0>(i), thrust::get<1>(i)}; + }); +} + template void bound_presolve_t::calc_and_set_updated_constraint_bounds(problem_t& pb) { diff --git a/cpp/src/mip/presolve/bounds_presolve.cuh b/cpp/src/mip/presolve/bounds_presolve.cuh index 84853a7812..0cb818e8ca 100644 --- a/cpp/src/mip/presolve/bounds_presolve.cuh +++ b/cpp/src/mip/presolve/bounds_presolve.cuh @@ -50,9 +50,7 @@ class bound_presolve_t { // when we need to accept a vector, we can use input_lb version termination_criterion_t solve(problem_t& pb, f_t lb, f_t ub, i_t var_idx); - termination_criterion_t solve(problem_t& pb, - raft::device_span input_lb = {}, - raft::device_span input_ub = {}); + termination_criterion_t solve(problem_t& pb); termination_criterion_t solve(problem_t& pb, const std::vector>& var_probe_val_pairs, @@ -62,6 +60,8 @@ class bound_presolve_t { void calculate_activity_on_problem_bounds(problem_t& pb); bool calculate_bounds_update(problem_t& pb); void set_updated_bounds(problem_t& pb); + void set_updated_bounds(const raft::handle_t* handle_ptr, + raft::device_span::type> output_bounds); void set_updated_bounds(const raft::handle_t* handle_ptr, raft::device_span output_lb, raft::device_span output_ub); diff --git a/cpp/src/mip/presolve/conditional_bound_strengthening.cu b/cpp/src/mip/presolve/conditional_bound_strengthening.cu index 2723a9b244..4eb8417b73 100644 --- a/cpp/src/mip/presolve/conditional_bound_strengthening.cu +++ b/cpp/src/mip/presolve/conditional_bound_strengthening.cu @@ -497,10 +497,10 @@ __global__ void update_constraint_bounds_kernel(typename problem_t::vi raft::device_span lock_per_constraint) { auto constraint_pair = constraint_pairs[blockIdx.x]; - int constr_i = constraint_pair.x; + int constr_i = get_lower(constraint_pair); if (constr_i == -1) { return; } - int constr_j = constraint_pair.y; + int constr_j = get_upper(constraint_pair); // FIXME:: for now handle only the constraints that fit in shared i_t offset_j = pb.offsets[constr_j]; @@ -550,8 +550,9 @@ __global__ void update_constraint_bounds_kernel(typename problem_t::vi if (tid < n_variables_in_constraint) { i_t variable_j = pb.variables[offset_j + tid]; a[tid] = pb.coefficients[offset_j + tid]; - lb[tid] = pb.variable_lower_bounds[variable_j]; - ub[tid] = pb.variable_upper_bounds[variable_j]; + auto bounds = pb.variable_bounds[variable_j]; + lb[tid] = get_lower(bounds); + ub[tid] = get_upper(bounds); vtypes[tid] = pb.variable_types[variable_j]; c[tid] = 0.; @@ -575,8 +576,9 @@ __global__ void update_constraint_bounds_kernel(typename problem_t::vi if (jj < 0) { f_t coeff = pb.coefficients[offset_i + index]; - f_t li = pb.variable_lower_bounds[variable_i]; - f_t ui = pb.variable_upper_bounds[variable_i]; + auto bounds = pb.variable_bounds[variable_i]; + f_t li = get_lower(bounds); + f_t ui = get_upper(bounds); min_activity_if_not_participating += (coeff > 0. ? coeff * li : coeff * ui); max_activity_if_not_participating += (coeff > 0. ? coeff * ui : coeff * li); } diff --git a/cpp/src/mip/presolve/lb_probing_cache.cu b/cpp/src/mip/presolve/lb_probing_cache.cu index 598a4c6bce..00c28d3450 100644 --- a/cpp/src/mip/presolve/lb_probing_cache.cu +++ b/cpp/src/mip/presolve/lb_probing_cache.cu @@ -201,9 +201,9 @@ __global__ void compute_min_slack_per_var(typename problem_t::view_t p if (std::signbit(a) != std::signbit(first_coeff)) { different_coeff = true; } auto cnst_idx = pb.reverse_constraints[var_offset + i]; auto cnstr_slack = cnst_slack[cnst_idx]; - auto delta_min_act = cnstr_slack.x + ((a < 0) ? a * ub : a * lb); + auto delta_min_act = get_lower(cnstr_slack) + ((a < 0) ? a * ub : a * lb); th_var_unit_slack = min(th_var_unit_slack, (delta_min_act / a)); - auto delta_max_act = cnstr_slack.y + ((a > 0) ? a * ub : a * lb); + auto delta_max_act = get_upper(cnstr_slack) + ((a > 0) ? a * ub : a * lb); th_var_unit_slack = min(th_var_unit_slack, (delta_max_act / a)); } __shared__ f_t shmem[raft::WarpSize]; @@ -232,7 +232,7 @@ __global__ void compute_min_slack_per_var(typename problem_t::view_t p th_max_excess = max(th_max_excess, excess); th_n_of_excess++; } - excess = max(0., cnstr_slack.y + diff); + excess = max(0., get_upper(cnstr_slack) + diff); if (excess > 0) { th_max_excess = max(th_max_excess, excess); th_n_of_excess++; diff --git a/cpp/src/mip/presolve/load_balanced_partition_helpers.cuh b/cpp/src/mip/presolve/load_balanced_partition_helpers.cuh index 55c18a902e..cffc5debb6 100644 --- a/cpp/src/mip/presolve/load_balanced_partition_helpers.cuh +++ b/cpp/src/mip/presolve/load_balanced_partition_helpers.cuh @@ -27,45 +27,6 @@ namespace cuopt::linear_programming::detail { -template -struct type_2 { - using type = void; -}; - -template <> -struct type_2 { - using type = int2; -}; - -template <> -struct type_2 { - using type = float2; -}; - -template <> -struct type_2 { - using type = double2; -}; - -template -raft::device_span::type> make_span_2(rmm::device_uvector& container) -{ - // TODO : ceildiv or throw assert - using T2 = typename type_2::type; - return raft::device_span(reinterpret_cast(container.data()), - sizeof(T) * container.size() / sizeof(T2)); -} - -template -raft::device_span::type> make_span_2( - rmm::device_uvector const& container) -{ - // TODO : ceildiv or throw assert - using T2 = typename type_2::type; - return raft::device_span(reinterpret_cast(container.data()), - sizeof(T) * container.size() / sizeof(T2)); -} - template constexpr int BitsPWrd = sizeof(degree_t) * 8; diff --git a/cpp/src/mip/presolve/multi_probe.cu b/cpp/src/mip/presolve/multi_probe.cu index 699a5f1dd4..ccd3f19511 100644 --- a/cpp/src/mip/presolve/multi_probe.cu +++ b/cpp/src/mip/presolve/multi_probe.cu @@ -331,36 +331,47 @@ void multi_probe_t::update_device_bounds(const raft::handle_t* handle_ } 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) +void multi_probe_t::update_host_bounds( + const raft::handle_t* handle_ptr, + const raft::device_span::type> variable_bounds) { - 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()); + cuopt_assert(variable_bounds.size() == host_lb.size(), "size of variable lower bound mismatch"); + cuopt_assert(variable_bounds.size() == host_ub.size(), "size of variable upper bound mismatch"); + + rmm::device_uvector var_lb(variable_bounds.size(), handle_ptr->get_stream()); + rmm::device_uvector var_ub(variable_bounds.size(), handle_ptr->get_stream()); + thrust::transform( + handle_ptr->get_thrust_policy(), + variable_bounds.begin(), + variable_bounds.end(), + thrust::make_zip_iterator(thrust::make_tuple(var_lb.begin(), var_ub.begin())), + [] __device__(auto i) { return thrust::make_tuple(get_lower(i), get_upper(i)); }); + raft::copy(host_lb.data(), var_lb.data(), var_lb.size(), handle_ptr->get_stream()); + raft::copy(host_ub.data(), var_ub.data(), var_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) { - cuopt_assert(upd_0.lb.size() == pb.variable_lower_bounds.size(), + cuopt_assert(upd_0.lb.size() == pb.variable_bounds.size(), "size of variable lower bound mismatch"); - raft::copy( - upd_0.lb.data(), pb.variable_lower_bounds.data(), upd_0.lb.size(), handle_ptr->get_stream()); - cuopt_assert(upd_1.lb.size() == pb.variable_lower_bounds.size(), + cuopt_assert(upd_1.lb.size() == pb.variable_bounds.size(), "size of variable lower bound mismatch"); - raft::copy( - upd_1.lb.data(), pb.variable_lower_bounds.data(), upd_1.lb.size(), handle_ptr->get_stream()); - cuopt_assert(upd_0.ub.size() == pb.variable_upper_bounds.size(), + cuopt_assert(upd_0.ub.size() == pb.variable_bounds.size(), "size of variable upper bound mismatch"); - raft::copy( - upd_0.ub.data(), pb.variable_upper_bounds.data(), upd_0.ub.size(), handle_ptr->get_stream()); - cuopt_assert(upd_1.ub.size() == pb.variable_upper_bounds.size(), + cuopt_assert(upd_1.ub.size() == pb.variable_bounds.size(), "size of variable upper bound mismatch"); - raft::copy( - upd_1.ub.data(), pb.variable_upper_bounds.data(), upd_1.ub.size(), handle_ptr->get_stream()); + + thrust::transform( + handle_ptr->get_thrust_policy(), + pb.variable_bounds.begin(), + pb.variable_bounds.end(), + thrust::make_zip_iterator( + thrust::make_tuple(upd_0.lb.begin(), upd_0.ub.begin(), upd_1.lb.begin(), upd_1.ub.begin())), + [] __device__(auto i) { + return thrust::make_tuple(get_lower(i), get_upper(i), get_lower(i), get_upper(i)); + }); } template @@ -410,6 +421,26 @@ void multi_probe_t::set_updated_bounds(const raft::handle_t* handle_pt raft::copy(output_lb.data(), lb.data(), lb.size(), handle_ptr->get_stream()); } +template +void multi_probe_t::set_updated_bounds( + const raft::handle_t* handle_ptr, + raft::device_span::type> output_bounds, + i_t select_update) +{ + auto& lb = select_update ? upd_1.lb : upd_0.lb; + auto& ub = select_update ? upd_1.ub : upd_0.ub; + + cuopt_assert(ub.size() == output_bounds.size(), "size of variable upper bound mismatch"); + cuopt_assert(lb.size() == output_bounds.size(), "size of variable lower bound mismatch"); + thrust::transform(handle_ptr->get_thrust_policy(), + thrust::make_zip_iterator(thrust::make_tuple(lb.begin(), ub.begin())), + thrust::make_zip_iterator(thrust::make_tuple(lb.end(), ub.end())), + output_bounds.begin(), + [] __device__(auto i) { + return typename type_2::type{thrust::get<0>(i), thrust::get<1>(i)}; + }); +} + template void multi_probe_t::constraint_stats(problem_t& pb, const raft::handle_t* handle_ptr) @@ -454,10 +485,7 @@ void multi_probe_t::set_updated_bounds(problem_t& pb, i_t select_update, const raft::handle_t* handle_ptr) { - set_updated_bounds(handle_ptr, - make_span(pb.variable_lower_bounds), - make_span(pb.variable_upper_bounds), - select_update); + set_updated_bounds(handle_ptr, make_span(pb.variable_bounds), select_update); } #if MIP_INSTANTIATE_FLOAT diff --git a/cpp/src/mip/presolve/multi_probe.cuh b/cpp/src/mip/presolve/multi_probe.cuh index 609702e909..db0b59f4ba 100644 --- a/cpp/src/mip/presolve/multi_probe.cuh +++ b/cpp/src/mip/presolve/multi_probe.cuh @@ -55,6 +55,9 @@ class multi_probe_t { void set_updated_bounds(problem_t& pb, i_t select_update, const raft::handle_t* handle_ptr); + void set_updated_bounds(const raft::handle_t* handle_ptr, + raft::device_span::type> output_bounds, + i_t select_update); void set_updated_bounds(const raft::handle_t* handle_ptr, raft::device_span output_lb, raft::device_span output_ub, @@ -72,8 +75,7 @@ class multi_probe_t { 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); + const raft::device_span::type> variable_bounds); void update_device_bounds(const raft::handle_t* handle_ptr); mip_solver_context_t& context; bounds_update_data_t upd_0; diff --git a/cpp/src/mip/presolve/probing_cache.cu b/cpp/src/mip/presolve/probing_cache.cu index 36140c5b31..7389ec206f 100644 --- a/cpp/src/mip/presolve/probing_cache.cu +++ b/cpp/src/mip/presolve/probing_cache.cu @@ -155,12 +155,11 @@ bool probing_cache_t::contains(problem_t& problem, i_t var_i return probing_cache.count(problem.original_ids[var_id]) > 0; } -template +template void inline insert_current_probing_to_cache(i_t var_idx, const val_interval_t& probe_val, bound_presolve_t& bound_presolve, - const std::vector& original_lb, - const std::vector& original_ub, + const std::vector& original_bounds, const std::vector& modified_lb, const std::vector& modified_ub, const std::vector& h_integer_indices, @@ -171,15 +170,16 @@ void inline insert_current_probing_to_cache(i_t var_idx, cache_entry_t cache_item; cache_item.val_interval = probe_val; for (auto impacted_var_idx : h_integer_indices) { - if (original_lb[impacted_var_idx] != modified_lb[impacted_var_idx] || - original_ub[impacted_var_idx] != modified_ub[impacted_var_idx]) { + auto original_var_bounds = original_bounds[impacted_var_idx]; + if (get_lower(original_var_bounds) != modified_lb[impacted_var_idx] || + get_upper(original_var_bounds) != modified_ub[impacted_var_idx]) { if (integer_equal( modified_lb[impacted_var_idx], modified_ub[impacted_var_idx], int_tol)) { ++n_implied_singletons; } - cuopt_assert(modified_lb[impacted_var_idx] >= original_lb[impacted_var_idx], + cuopt_assert(modified_lb[impacted_var_idx] >= get_lower(original_var_bounds), "Lower bound must be greater than or equal to original lower bound"); - cuopt_assert(modified_ub[impacted_var_idx] <= original_ub[impacted_var_idx], + cuopt_assert(modified_ub[impacted_var_idx] <= get_upper(original_var_bounds), "Upper bound must be less than or equal to original upper bound"); cached_bound_t new_bound{modified_lb[impacted_var_idx], modified_ub[impacted_var_idx]}; cache_item.var_to_cached_bound_map.insert({impacted_var_idx, new_bound}); @@ -210,8 +210,9 @@ __global__ void compute_min_slack_per_var(typename problem_t::view_t p i_t var_offset = pb.reverse_offsets[var_idx]; i_t var_degree = pb.reverse_offsets[var_idx + 1] - var_offset; f_t th_var_unit_slack = std::numeric_limits::max(); - f_t lb = pb.variable_lower_bounds[var_idx]; - f_t ub = pb.variable_upper_bounds[var_idx]; + auto var_bounds = pb.variable_bounds[var_idx]; + f_t lb = get_lower(var_bounds); + f_t ub = get_upper(var_bounds); f_t first_coeff = pb.reverse_coefficients[var_offset]; bool different_coeff = false; for (i_t i = threadIdx.x; i < var_degree; i += blockDim.x) { @@ -360,13 +361,12 @@ inline std::vector compute_prioritized_integer_indices( return h_priority_indices; } -template +template void compute_cache_for_var(i_t var_idx, bound_presolve_t& bound_presolve, problem_t& problem, multi_probe_t& multi_probe_presolve, - const std::vector& h_var_lower_bounds, - const std::vector& h_var_upper_bounds, + const std::vector& h_var_bounds, const std::vector& h_integer_indices, std::atomic& n_of_implied_singletons, std::atomic& n_of_cached_probings, @@ -375,11 +375,12 @@ void compute_cache_for_var(i_t var_idx, RAFT_CUDA_TRY(cudaSetDevice(device_id)); // test if we need per thread handle raft::handle_t handle{}; - std::vector h_improved_lower_bounds(h_var_lower_bounds.size()); - std::vector h_improved_upper_bounds(h_var_upper_bounds.size()); + std::vector h_improved_lower_bounds(h_var_bounds.size()); + std::vector h_improved_upper_bounds(h_var_bounds.size()); std::pair, val_interval_t> probe_vals; - f_t lb = h_var_lower_bounds[var_idx]; - f_t ub = h_var_upper_bounds[var_idx]; + auto bounds = h_var_bounds[var_idx]; + f_t lb = get_lower(bounds); + f_t ub = get_upper(bounds); for (i_t i = 0; i < 2; ++i) { auto& probe_val = i == 0 ? probe_vals.first : probe_vals.second; // if binary, probe both values @@ -451,8 +452,7 @@ void compute_cache_for_var(i_t var_idx, insert_current_probing_to_cache(var_idx, probe_val, bound_presolve, - h_var_lower_bounds, - h_var_upper_bounds, + h_var_bounds, h_improved_lower_bounds, h_improved_upper_bounds, h_integer_indices, @@ -470,9 +470,8 @@ void compute_probing_cache(bound_presolve_t& bound_presolve, // we dont want to compute the probing cache for all variables for time and computation resources auto priority_indices = compute_prioritized_integer_indices(bound_presolve, problem); CUOPT_LOG_DEBUG("Computing probing cache"); - auto h_integer_indices = host_copy(problem.integer_indices); - const auto h_var_upper_bounds = host_copy(problem.variable_upper_bounds); - const auto h_var_lower_bounds = host_copy(problem.variable_lower_bounds); + auto h_integer_indices = host_copy(problem.integer_indices); + const auto h_var_bounds = host_copy(problem.variable_bounds); // TODO adjust the iteration limit depending on the total time limit and time it takes for single // var bound_presolve.settings.iteration_limit = 50; @@ -512,8 +511,7 @@ void compute_probing_cache(bound_presolve_t& bound_presolve, bound_presolve, problem, multi_probe_presolve, - h_var_lower_bounds, - h_var_upper_bounds, + h_var_bounds, h_integer_indices, n_of_implied_singletons, n_of_cached_probings, diff --git a/cpp/src/mip/presolve/probing_cache.cuh b/cpp/src/mip/presolve/probing_cache.cuh index 2ba5010c6b..755c18b0bb 100644 --- a/cpp/src/mip/presolve/probing_cache.cuh +++ b/cpp/src/mip/presolve/probing_cache.cuh @@ -17,8 +17,6 @@ #pragma once -#include -#include #include "bounds_presolve.cuh" #include @@ -30,12 +28,6 @@ namespace cuopt::linear_programming::detail { template class bound_presolve_t; -template -class load_balanced_bounds_presolve_t; - -template -class load_balanced_problem_t; - /* Probing cache is a set of implied bounds when we set a variable to some value. We keep two sets of changed bounds for each interval: @@ -137,9 +129,4 @@ void compute_probing_cache(bound_presolve_t& bound_presolve, problem_t& problem, timer_t timer); -template -void compute_probing_cache(load_balanced_bounds_presolve_t& bound_presolve, - load_balanced_problem_t& problem, - timer_t timer); - } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/trivial_presolve.cuh b/cpp/src/mip/presolve/trivial_presolve.cuh index f4940a62f6..803bfca800 100644 --- a/cpp/src/mip/presolve/trivial_presolve.cuh +++ b/cpp/src/mip/presolve/trivial_presolve.cuh @@ -77,16 +77,11 @@ void cleanup_vectors(problem_t& pb, handle_ptr->get_stream()); handle_ptr->sync_stream(); - auto lb_iter = thrust::remove_if(handle_ptr->get_thrust_policy(), - pb.variable_lower_bounds.begin(), - pb.variable_lower_bounds.end(), - var_map.begin(), - is_zero_t{}); - auto ub_iter = thrust::remove_if(handle_ptr->get_thrust_policy(), - pb.variable_upper_bounds.begin(), - pb.variable_upper_bounds.end(), - var_map.begin(), - is_zero_t{}); + auto bnd_iter = thrust::remove_if(handle_ptr->get_thrust_policy(), + pb.variable_bounds.begin(), + pb.variable_bounds.end(), + var_map.begin(), + is_zero_t{}); auto type_iter = thrust::remove_if(handle_ptr->get_thrust_policy(), pb.variable_types.begin(), pb.variable_types.end(), @@ -102,10 +97,7 @@ void cleanup_vectors(problem_t& pb, pb.objective_coefficients.end(), var_map.begin(), is_zero_t{}); - pb.variable_lower_bounds.resize(lb_iter - pb.variable_lower_bounds.begin(), - handle_ptr->get_stream()); - pb.variable_upper_bounds.resize(ub_iter - pb.variable_upper_bounds.begin(), - handle_ptr->get_stream()); + pb.variable_bounds.resize(bnd_iter - pb.variable_bounds.begin(), handle_ptr->get_stream()); pb.variable_types.resize(type_iter - pb.variable_types.begin(), handle_ptr->get_stream()); pb.is_binary_variable.resize(binary_iter - pb.is_binary_variable.begin(), handle_ptr->get_stream()); @@ -117,6 +109,7 @@ void cleanup_vectors(problem_t& pb, template void update_from_csr(problem_t& pb) { + using f_t2 = typename type_2::type; auto handle_ptr = pb.handle_ptr; rmm::device_uvector cnst(pb.coefficients.size(), handle_ptr->get_stream()); thrust::uninitialized_fill(handle_ptr->get_thrust_policy(), cnst.begin(), cnst.end(), 0); @@ -145,9 +138,8 @@ void update_from_csr(problem_t& pb) thrust::stable_partition(handle_ptr->get_thrust_policy(), coo_begin, coo_begin + cnst.size(), - is_variable_free_t{pb.tolerances.integrality_tolerance, - make_span(pb.variable_lower_bounds), - make_span(pb.variable_upper_bounds)}); + is_variable_free_t{pb.tolerances.integrality_tolerance, + make_span(pb.variable_bounds)}); RAFT_CHECK_CUDA(handle_ptr->get_stream()); nnz_edge_count = partition_iter - coo_begin; } @@ -180,12 +172,11 @@ void update_from_csr(problem_t& pb) handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), thrust::make_counting_iterator(pb.n_variables), - assign_fixed_var_t{make_span(var_map), - make_span(pb.variable_lower_bounds), - make_span(pb.variable_upper_bounds), - make_span(pb.objective_coefficients), - make_span(pb.presolve_data.variable_mapping), - make_span(pb.presolve_data.fixed_var_assignment)}); + assign_fixed_var_t{make_span(var_map), + make_span(pb.variable_bounds), + make_span(pb.objective_coefficients), + make_span(pb.presolve_data.variable_mapping), + make_span(pb.presolve_data.fixed_var_assignment)}); auto used_iter = thrust::stable_partition(handle_ptr->get_thrust_policy(), pb.presolve_data.variable_mapping.begin(), pb.presolve_data.variable_mapping.end(), @@ -203,11 +194,10 @@ void update_from_csr(problem_t& pb) handle_ptr->get_stream()); rmm::device_uvector unused_coo_cnst_bound_updates(cnst.size() - nnz_edge_count, handle_ptr->get_stream()); - elem_multi_t mul{make_span(pb.coefficients), - make_span(pb.variables), - make_span(pb.objective_coefficients), - make_span(pb.variable_lower_bounds), - make_span(pb.variable_upper_bounds)}; + elem_multi_t mul{make_span(pb.coefficients), + make_span(pb.variables), + make_span(pb.objective_coefficients), + make_span(pb.variable_bounds)}; auto iter = thrust::reduce_by_key( handle_ptr->get_thrust_policy(), @@ -233,16 +223,14 @@ void update_from_csr(problem_t& pb) } // update objective_offset - pb.presolve_data.objective_offset += - thrust::transform_reduce(handle_ptr->get_thrust_policy(), - thrust::counting_iterator(0), - thrust::counting_iterator(pb.n_variables), - unused_var_obj_offset_t{make_span(var_map), - make_span(pb.objective_coefficients), - make_span(pb.variable_lower_bounds), - make_span(pb.variable_upper_bounds)}, - 0., - thrust::plus{}); + pb.presolve_data.objective_offset += thrust::transform_reduce( + handle_ptr->get_thrust_policy(), + thrust::counting_iterator(0), + thrust::counting_iterator(pb.n_variables), + unused_var_obj_offset_t{ + make_span(var_map), make_span(pb.objective_coefficients), make_span(pb.variable_bounds)}, + 0., + thrust::plus{}); RAFT_CHECK_CUDA(handle_ptr->get_stream()); // create renumbering maps diff --git a/cpp/src/mip/presolve/trivial_presolve_helpers.cuh b/cpp/src/mip/presolve/trivial_presolve_helpers.cuh index b7e3f4b464..cf7e5064d2 100644 --- a/cpp/src/mip/presolve/trivial_presolve_helpers.cuh +++ b/cpp/src/mip/presolve/trivial_presolve_helpers.cuh @@ -30,40 +30,34 @@ struct non_zero_degree_t { __device__ i_t operator()(i_t i) { return offsets[i] != offsets[i + 1]; } }; -template +template struct is_variable_free_t { f_t tol; - raft::device_span lb; - raft::device_span ub; - is_variable_free_t(f_t tol_, raft::device_span lb_, raft::device_span ub_) - : tol(tol_), lb(lb_), ub(ub_) - { - } + raft::device_span bnd; + is_variable_free_t(f_t tol_, raft::device_span bnd_) : tol(tol_), bnd(bnd_) {} template __device__ bool operator()(tuple_t edge) { - auto var = thrust::get<2>(edge); - return abs(ub[var] - lb[var]) > tol; + auto var = thrust::get<2>(edge); + auto bounds = bnd[var]; + return abs(get_upper(bounds) - get_lower(bounds)) > tol; } }; -template +template struct assign_fixed_var_t { raft::device_span is_var_used; - raft::device_span variable_lower_bounds; - raft::device_span variable_upper_bounds; + raft::device_span variable_bounds; raft::device_span objective_coefficients; raft::device_span variable_mapping; raft::device_span fixed_assignment; assign_fixed_var_t(raft::device_span is_var_used_, - raft::device_span variable_lower_bounds_, - raft::device_span variable_upper_bounds_, + raft::device_span variable_bounds_, raft::device_span objective_coefficients_, raft::device_span variable_mapping_, raft::device_span fixed_assignment_) : is_var_used(is_var_used_), - variable_lower_bounds(variable_lower_bounds_), - variable_upper_bounds(variable_upper_bounds_), + variable_bounds(variable_bounds_), objective_coefficients(objective_coefficients_), variable_mapping(variable_mapping_), fixed_assignment(fixed_assignment_) @@ -74,39 +68,38 @@ struct assign_fixed_var_t { { if (!is_var_used[i]) { auto orig_v_idx = variable_mapping[i]; + auto bounds = variable_bounds[i]; fixed_assignment[orig_v_idx] = - (objective_coefficients[i] > 0) ? variable_lower_bounds[i] : variable_upper_bounds[i]; + (objective_coefficients[i] > 0) ? get_lower(bounds) : get_upper(bounds); } } }; -template +template struct elem_multi_t { raft::device_span coefficients; raft::device_span variables; raft::device_span obj_coefficients; - raft::device_span variable_lower_bounds; - raft::device_span variable_upper_bounds; + raft::device_span variable_bounds; elem_multi_t(raft::device_span coefficients_, raft::device_span variables_, raft::device_span obj_coefficients_, - raft::device_span variable_lower_bounds_, - raft::device_span variable_upper_bounds_) + raft::device_span variable_bounds_) : coefficients(coefficients_), variables(variables_), obj_coefficients(obj_coefficients_), - variable_lower_bounds(variable_lower_bounds_), - variable_upper_bounds(variable_upper_bounds_) + variable_bounds(variable_bounds_) { } __device__ f_t operator()(i_t i) const { - auto var = variables[i]; + auto var = variables[i]; + auto bounds = variable_bounds[var]; if (obj_coefficients[var] > 0) { - return variable_lower_bounds[var] * coefficients[i]; + return get_lower(bounds) * coefficients[i]; } else { - return variable_upper_bounds[var] * coefficients[i]; + return get_upper(bounds) * coefficients[i]; } } }; @@ -136,18 +129,16 @@ struct update_constraint_bounds_t { } }; -template +template struct unused_var_obj_offset_t { raft::device_span var_map; raft::device_span objective_coefficients; - raft::device_span lb; - raft::device_span ub; + raft::device_span bnd; unused_var_obj_offset_t(raft::device_span var_map_, raft::device_span objective_coefficients_, - raft::device_span lb_, - raft::device_span ub_) - : var_map(var_map_), objective_coefficients(objective_coefficients_), lb(lb_), ub(ub_) + raft::device_span bnd_) + : var_map(var_map_), objective_coefficients(objective_coefficients_), bnd(bnd_) { } @@ -156,7 +147,8 @@ struct unused_var_obj_offset_t { auto obj_coeff = objective_coefficients[i]; // in case both bounds are infinite if (obj_coeff == 0.) return 0.; - auto obj_off = (obj_coeff > 0) ? obj_coeff * lb[i] : obj_coeff * ub[i]; + auto bounds = bnd[i]; + auto obj_off = (obj_coeff > 0) ? obj_coeff * get_lower(bounds) : obj_coeff * get_upper(bounds); return var_map[i] ? 0. : obj_off; } }; diff --git a/cpp/src/mip/problem/host_helper.cuh b/cpp/src/mip/problem/host_helper.cuh index 4023a32d27..9136e6899e 100644 --- a/cpp/src/mip/problem/host_helper.cuh +++ b/cpp/src/mip/problem/host_helper.cuh @@ -17,6 +17,7 @@ #pragma once +#include #include #include @@ -52,22 +53,21 @@ struct constraints_delta_t { template struct variables_delta_t { + using f_t2 = typename type_2::type; std::vector objective_coefficients; - std::vector lower_bounds; - std::vector upper_bounds; + std::vector variable_bounds; std::vector variable_types; std::vector is_binary_variable; i_t n_vars; - i_t size() const { return lower_bounds.size(); } + i_t size() const { return variable_bounds.size(); } // returns the added variable id i_t add_variable(f_t lower_bound, f_t upper_bound, f_t obj_weight, var_t var_type) { cuopt_assert(lower_bound >= 0, "Variable bounds must be non-negative!"); - lower_bounds.push_back(lower_bound); - upper_bounds.push_back(upper_bound); + variable_bounds.push_back(f_t2{lower_bound, upper_bound}); objective_coefficients.push_back(obj_weight); variable_types.push_back(var_type); is_binary_variable.push_back(0); diff --git a/cpp/src/mip/problem/problem.cu b/cpp/src/mip/problem/problem.cu index 1a5f76b038..37064b5c63 100644 --- a/cpp/src/mip/problem/problem.cu +++ b/cpp/src/mip/problem/problem.cu @@ -70,6 +70,8 @@ void problem_t::op_problem_cstr_body(const optimization_problem_tget_problem_category() != problem_category_t::LP; if (is_mip) { variable_types = @@ -121,10 +123,7 @@ problem_t::problem_t( offsets(problem_.get_constraint_matrix_offsets(), problem_.get_handle_ptr()->get_stream()), objective_coefficients(problem_.get_objective_coefficients(), problem_.get_handle_ptr()->get_stream()), - variable_lower_bounds(problem_.get_variable_lower_bounds(), - problem_.get_handle_ptr()->get_stream()), - variable_upper_bounds(problem_.get_variable_upper_bounds(), - problem_.get_handle_ptr()->get_stream()), + variable_bounds(0, problem_.get_handle_ptr()->get_stream()), constraint_lower_bounds(problem_.get_constraint_lower_bounds(), problem_.get_handle_ptr()->get_stream()), constraint_upper_bounds(problem_.get_constraint_upper_bounds(), @@ -173,8 +172,7 @@ problem_t::problem_t(const problem_t& problem_) variables(problem_.variables, handle_ptr->get_stream()), offsets(problem_.offsets, handle_ptr->get_stream()), objective_coefficients(problem_.objective_coefficients, handle_ptr->get_stream()), - variable_lower_bounds(problem_.variable_lower_bounds, handle_ptr->get_stream()), - variable_upper_bounds(problem_.variable_upper_bounds, handle_ptr->get_stream()), + variable_bounds(problem_.variable_bounds, handle_ptr->get_stream()), constraint_lower_bounds(problem_.constraint_lower_bounds, handle_ptr->get_stream()), constraint_upper_bounds(problem_.constraint_upper_bounds, handle_ptr->get_stream()), combined_bounds(problem_.combined_bounds, handle_ptr->get_stream()), @@ -247,16 +245,10 @@ problem_t::problem_t(const problem_t& problem_, bool no_deep ? rmm::device_uvector(problem_.objective_coefficients, handle_ptr->get_stream()) : rmm::device_uvector(problem_.objective_coefficients.size(), handle_ptr->get_stream())), - variable_lower_bounds( - (!no_deep_copy) - ? rmm::device_uvector(problem_.variable_lower_bounds, handle_ptr->get_stream()) - : rmm::device_uvector(problem_.variable_lower_bounds.size(), - handle_ptr->get_stream())), - variable_upper_bounds( + variable_bounds( (!no_deep_copy) - ? rmm::device_uvector(problem_.variable_upper_bounds, handle_ptr->get_stream()) - : rmm::device_uvector(problem_.variable_upper_bounds.size(), - handle_ptr->get_stream())), + ? rmm::device_uvector(problem_.variable_bounds, handle_ptr->get_stream()) + : rmm::device_uvector(problem_.variable_bounds.size(), handle_ptr->get_stream())), constraint_lower_bounds( (!no_deep_copy) ? rmm::device_uvector(problem_.constraint_lower_bounds, handle_ptr->get_stream()) @@ -389,16 +381,12 @@ void problem_t::check_problem_representation(bool check_transposed, } // Check variable bounds are set and with the correct size - if (!empty) { - cuopt_assert(!variable_lower_bounds.is_empty() && !variable_upper_bounds.is_empty(), - "Variable lower bounds and variable upper bounds must be set."); - } - cuopt_assert(variable_lower_bounds.size() == objective_coefficients.size(), + if (!empty) { cuopt_assert(!variable_bounds.is_empty(), "Variable bounds must be set."); } + cuopt_assert(variable_bounds.size() == objective_coefficients.size(), "Sizes for vectors related to the variables are not the same."); - cuopt_assert(variable_upper_bounds.size() == objective_coefficients.size(), - "Sizes for vectors related to the variables are not the same"); - cuopt_assert(variable_upper_bounds.size() == (std::size_t)n_variables, + cuopt_assert(variable_bounds.size() == (std::size_t)n_variables, "Sizes for vectors related to the variables are not the same."); + cuopt_assert(variable_types.size() == (std::size_t)n_variables, "Sizes for vectors related to the variables are not the same."); // Check constraints bounds sizes @@ -418,16 +406,15 @@ void problem_t::check_problem_representation(bool check_transposed, "Sizes for vectors related to the constraints are not the same."); // Check the validity of bounds - cuopt_expects( - thrust::all_of(handle_ptr->get_thrust_policy(), - thrust::make_counting_iterator(0), - thrust::make_counting_iterator(n_variables), - [variable_lower_bounds = variable_lower_bounds.data(), - variable_upper_bounds = variable_upper_bounds.data()] __device__(i_t idx) { - return variable_lower_bounds[idx] <= variable_upper_bounds[idx]; - }), - error_type_t::ValidationError, - "Variable bounds are invalid"); + cuopt_expects(thrust::all_of(handle_ptr->get_thrust_policy(), + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(n_variables), + [vars_bnd = make_span(variable_bounds)] __device__(i_t idx) { + auto bounds = vars_bnd[idx]; + return get_lower(bounds) <= get_upper(bounds); + }), + error_type_t::ValidationError, + "Variable bounds are invalid"); cuopt_expects( thrust::all_of(handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), @@ -574,23 +561,23 @@ void problem_t::check_problem_representation(bool check_transposed, return true; }), "Some variables aren't referenced in the appropriate indice tables"); - cuopt_assert( - thrust::all_of( - handle_ptr->get_thrust_policy(), - thrust::make_zip_iterator(thrust::make_counting_iterator(0), - is_binary_variable.cbegin()), - thrust::make_zip_iterator(thrust::make_counting_iterator(is_binary_variable.size()), + cuopt_assert(thrust::all_of(handle_ptr->get_thrust_policy(), + thrust::make_zip_iterator(thrust::make_counting_iterator(0), + is_binary_variable.cbegin()), + thrust::make_zip_iterator( + thrust::make_counting_iterator(is_binary_variable.size()), is_binary_variable.cend()), - [types = variable_types.data(), - lb = variable_lower_bounds.data(), - ub = variable_upper_bounds.data(), - v = view()] __device__(const thrust::tuple tuple) { - i_t idx = thrust::get<0>(tuple); - i_t pred = thrust::get<1>(tuple); - return pred == (types[idx] != var_t::CONTINUOUS && v.integer_equal(lb[idx], 0.) && - v.integer_equal(ub[idx], 1.)); - }), - "The binary variable table is incorrect."); + [types = variable_types.data(), + vars_bnd = make_span(variable_bounds), + v = view()] __device__(const thrust::tuple tuple) { + i_t idx = thrust::get<0>(tuple); + i_t pred = thrust::get<1>(tuple); + auto bounds = vars_bnd[idx]; + return pred == (types[idx] != var_t::CONTINUOUS && + v.integer_equal(get_lower(bounds), 0.) && + v.integer_equal(get_upper(bounds), 1.)); + }), + "The binary variable table is incorrect."); if (!empty) { cuopt_assert(is_binary_pb == (n_variables == thrust::count(handle_ptr->get_thrust_policy(), is_binary_variable.begin(), @@ -768,9 +755,10 @@ void problem_t::compute_binary_var_table() is_binary_variable.begin(), is_binary_variable.end(), [pb_view] __device__(i_t i) { + auto bounds = pb_view.variable_bounds[i]; return pb_view.variable_types[i] != var_t::CONTINUOUS && - (pb_view.integer_equal(pb_view.variable_lower_bounds[i], 0) && - pb_view.integer_equal(pb_view.variable_upper_bounds[i], 1)); + (pb_view.integer_equal(get_lower(bounds), 0) && + pb_view.integer_equal(get_upper(bounds), 1)); }); get_n_binary_variables(); @@ -919,10 +907,7 @@ typename problem_t::view_t problem_t::view() v.offsets = raft::device_span{offsets.data(), offsets.size()}; v.objective_coefficients = raft::device_span{objective_coefficients.data(), objective_coefficients.size()}; - v.variable_lower_bounds = - raft::device_span{variable_lower_bounds.data(), variable_lower_bounds.size()}; - v.variable_upper_bounds = - raft::device_span{variable_upper_bounds.data(), variable_upper_bounds.size()}; + v.variable_bounds = make_span(variable_bounds); v.constraint_lower_bounds = raft::device_span{constraint_lower_bounds.data(), constraint_lower_bounds.size()}; v.constraint_upper_bounds = @@ -945,8 +930,7 @@ typename problem_t::view_t problem_t::view() template void problem_t::resize_variables(size_t size) { - variable_lower_bounds.resize(size, handle_ptr->get_stream()); - variable_upper_bounds.resize(size, handle_ptr->get_stream()); + variable_bounds.resize(size, handle_ptr->get_stream()); variable_types.resize(size, handle_ptr->get_stream()); objective_coefficients.resize(size, handle_ptr->get_stream()); is_binary_variable.resize(size, handle_ptr->get_stream()); @@ -978,13 +962,9 @@ void problem_t::insert_variables(variables_delta_t& h_vars) CUOPT_LOG_DEBUG("problem added variable size %d prev %d", h_vars.size(), n_variables); // resize the variable arrays if it can't fit the variables resize_variables(n_variables + h_vars.size()); - raft::copy(variable_lower_bounds.data() + n_variables, - h_vars.lower_bounds.data(), - h_vars.lower_bounds.size(), - handle_ptr->get_stream()); - raft::copy(variable_upper_bounds.data() + n_variables, - h_vars.upper_bounds.data(), - h_vars.upper_bounds.size(), + raft::copy(variable_bounds.data() + n_variables, + h_vars.variable_bounds.data(), + h_vars.variable_bounds.size(), handle_ptr->get_stream()); raft::copy(variable_types.data() + n_variables, h_vars.variable_types.data(), @@ -1225,15 +1205,9 @@ void problem_t::remove_given_variables(problem_t& original_p thrust::gather(handle_ptr->get_thrust_policy(), variable_map.begin(), variable_map.end(), - original_problem.variable_lower_bounds.begin(), - variable_lower_bounds.begin()); - variable_lower_bounds.resize(variable_map.size(), handle_ptr->get_stream()); - thrust::gather(handle_ptr->get_thrust_policy(), - variable_map.begin(), - variable_map.end(), - original_problem.variable_upper_bounds.begin(), - variable_upper_bounds.begin()); - variable_upper_bounds.resize(variable_map.size(), handle_ptr->get_stream()); + original_problem.variable_bounds.begin(), + variable_bounds.begin()); + variable_bounds.resize(variable_map.size(), handle_ptr->get_stream()); thrust::gather(handle_ptr->get_thrust_policy(), variable_map.begin(), variable_map.end(), @@ -1351,20 +1325,20 @@ void standardize_bounds(std::vector>>& variable_ problem_t& pb) { auto handle_ptr = pb.handle_ptr; - auto h_var_lower_bounds = cuopt::host_copy(pb.variable_lower_bounds); - auto h_var_upper_bounds = cuopt::host_copy(pb.variable_upper_bounds); + auto h_var_bounds = cuopt::host_copy(pb.variable_bounds); auto h_objective_coefficients = cuopt::host_copy(pb.objective_coefficients); auto h_variable_types = cuopt::host_copy(pb.variable_types); handle_ptr->sync_stream(); - const i_t n_vars_originally = (i_t)h_var_lower_bounds.size(); + const i_t n_vars_originally = (i_t)h_var_bounds.size(); for (i_t i = 0; i < n_vars_originally; ++i) { // if variable has free bounds, replace it with two vars // but add only one var and use it in all constraints // TODO create one var for integrals and one var for continuous - if (h_var_lower_bounds[i] == -std::numeric_limits::infinity() && - h_var_upper_bounds[i] == std::numeric_limits::infinity()) { + auto h_var_bound = h_var_bounds[i]; + if (get_lower(h_var_bound) == -std::numeric_limits::infinity() && + get_upper(h_var_bound) == std::numeric_limits::infinity()) { // add new variable auto var_coeff_vec = variable_constraint_map[i]; // negate all values in vec @@ -1372,16 +1346,16 @@ void standardize_bounds(std::vector>>& variable_ coeff = -coeff; } - h_var_lower_bounds[i] = 0.; + h_var_bounds[i].x = 0.; pb.presolve_data.variable_offsets[i] = 0.; pb.presolve_data.additional_var_used[i] = true; pb.presolve_data.additional_var_id_per_var[i] = pb.n_variables; + using f_t2 = typename type_2::type; // new var data std::stable_sort(var_coeff_vec.begin(), var_coeff_vec.end()); variable_constraint_map.push_back(var_coeff_vec); - h_var_lower_bounds.push_back(0.); - h_var_upper_bounds.push_back(std::numeric_limits::infinity()); + h_var_bounds.push_back(f_t2{0., std::numeric_limits::infinity()}); pb.presolve_data.variable_offsets.push_back(0.); h_objective_coefficients.push_back(-h_objective_coefficients[i]); h_variable_types.push_back(h_variable_types[i]); @@ -1397,21 +1371,14 @@ void standardize_bounds(std::vector>>& variable_ // TODO add some tests // resize the device vectors is sizes are smaller - if (pb.variable_lower_bounds.size() < h_var_lower_bounds.size()) { - pb.variable_lower_bounds.resize(h_var_lower_bounds.size(), handle_ptr->get_stream()); - pb.variable_upper_bounds.resize(h_var_lower_bounds.size(), handle_ptr->get_stream()); + if (pb.variable_bounds.size() < h_var_bounds.size()) { + pb.variable_bounds.resize(h_var_bounds.size(), handle_ptr->get_stream()); pb.objective_coefficients.resize(h_objective_coefficients.size(), handle_ptr->get_stream()); pb.variable_types.resize(h_variable_types.size(), handle_ptr->get_stream()); } - raft::copy(pb.variable_lower_bounds.data(), - h_var_lower_bounds.data(), - h_var_lower_bounds.size(), - handle_ptr->get_stream()); - raft::copy(pb.variable_upper_bounds.data(), - h_var_upper_bounds.data(), - h_var_upper_bounds.size(), - handle_ptr->get_stream()); + raft::copy( + pb.variable_bounds.data(), h_var_bounds.data(), h_var_bounds.size(), handle_ptr->get_stream()); raft::copy(pb.objective_coefficients.data(), h_objective_coefficients.data(), h_objective_coefficients.size(), @@ -1539,9 +1506,9 @@ void problem_t::get_host_user_problem( } } user_problem.num_range_rows = user_problem.range_rows.size(); - user_problem.lower = cuopt::host_copy(variable_lower_bounds); - user_problem.upper = cuopt::host_copy(variable_upper_bounds); - user_problem.problem_name = original_problem_ptr->get_problem_name(); + std::tie(user_problem.lower, user_problem.upper) = + extract_host_bounds(variable_bounds, handle_ptr); + user_problem.problem_name = original_problem_ptr->get_problem_name(); if (static_cast(row_names.size()) == m) { user_problem.row_names.resize(m); for (int i = 0; i < m; ++i) { diff --git a/cpp/src/mip/problem/problem.cuh b/cpp/src/mip/problem/problem.cuh index 49103fc955..9d63e18579 100644 --- a/cpp/src/mip/problem/problem.cuh +++ b/cpp/src/mip/problem/problem.cuh @@ -136,8 +136,9 @@ class problem_t { DI bool check_variable_within_bounds(i_t v, f_t val) const { const f_t int_tol = tolerances.integrality_tolerance; + auto bounds = variable_bounds[v]; bool within_bounds = - val <= (variable_upper_bounds[v] + int_tol) && val >= (variable_lower_bounds[v] - int_tol); + val <= (get_upper(bounds) + int_tol) && val >= (get_lower(bounds) - int_tol); return within_bounds; } @@ -158,21 +159,21 @@ class problem_t { { cuopt_assert(var_t::INTEGER != variable_types[v], "Random value can only be called on continuous values"); - f_t lower_bound = variable_lower_bounds[v]; - f_t upper_bound = variable_upper_bounds[v]; + auto bounds = variable_bounds[v]; f_t val; - if (isfinite(lower_bound) && isfinite(upper_bound)) { - f_t diff = upper_bound - lower_bound; - val = diff * rng.next_float() + lower_bound; + if (isfinite(get_lower(bounds)) && isfinite(get_upper(bounds))) { + f_t diff = get_upper(bounds) - get_lower(bounds); + val = diff * rng.next_float() + get_lower(bounds); } else { - auto finite_bound = isfinite(lower_bound) ? lower_bound : upper_bound; + auto finite_bound = isfinite(get_lower(bounds)) ? get_lower(bounds) : get_upper(bounds); val = finite_bound; } - cuopt_assert(isfinite(lower_bound), "Value must be finite"); + cuopt_assert(isfinite(get_lower(bounds)), "Value must be finite"); return val; } + using f_t2 = typename type_2::type; typename mip_solver_settings_t::tolerances_t tolerances; i_t n_variables; i_t n_integer_vars; @@ -187,8 +188,7 @@ class problem_t { raft::device_span variables; raft::device_span offsets; raft::device_span objective_coefficients; - raft::device_span variable_lower_bounds; - raft::device_span variable_upper_bounds; + raft::device_span variable_bounds; raft::device_span constraint_lower_bounds; raft::device_span constraint_upper_bounds; raft::device_span variable_types; @@ -242,8 +242,8 @@ class problem_t { /** weights in the objective function */ rmm::device_uvector objective_coefficients; - rmm::device_uvector variable_lower_bounds; - rmm::device_uvector variable_upper_bounds; + using f_t2 = typename type_2::type; + rmm::device_uvector variable_bounds; rmm::device_uvector constraint_lower_bounds; rmm::device_uvector constraint_upper_bounds; /* biggest between cstr lower and upper */ diff --git a/cpp/src/mip/problem/problem_helpers.cuh b/cpp/src/mip/problem/problem_helpers.cuh index bdc4635c03..8451983a99 100644 --- a/cpp/src/mip/problem/problem_helpers.cuh +++ b/cpp/src/mip/problem/problem_helpers.cuh @@ -58,6 +58,36 @@ struct transform_bounds_functor { } }; +template +static void set_variable_bounds(detail::problem_t& op_problem) +{ + op_problem.variable_bounds.resize(op_problem.n_variables, op_problem.handle_ptr->get_stream()); + auto vars_bnd = make_span(op_problem.variable_bounds); + + auto orig_problem = op_problem.original_problem_ptr; + auto variable_lower_bounds = make_span(orig_problem->get_variable_lower_bounds()); + auto variable_upper_bounds = make_span(orig_problem->get_variable_upper_bounds()); + + bool default_variable_lb = (orig_problem->get_variable_lower_bounds().is_empty()); + bool default_variable_ub = (orig_problem->get_variable_upper_bounds().is_empty()); + + thrust::for_each(op_problem.handle_ptr->get_thrust_policy(), + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(op_problem.n_variables), + [vars_bnd, + variable_lower_bounds, + variable_upper_bounds, + default_variable_lb, + default_variable_ub] __device__(auto i) { + using f_t2 = typename type_2::type; + auto lb = f_t{0}; + auto ub = std::numeric_limits::infinity(); + if (!default_variable_lb) { lb = variable_lower_bounds[i]; } + if (!default_variable_ub) { ub = variable_upper_bounds[i]; } + vars_bnd[i] = f_t2{lb, ub}; + }); +} + template static void set_bounds_if_not_set(detail::problem_t& op_problem) { @@ -90,25 +120,7 @@ static void set_bounds_if_not_set(detail::problem_t& op_problem) transform_bounds_functor()); } - // If variable bound was not set, set it to default value - if (op_problem.variable_lower_bounds.is_empty() && - !op_problem.objective_coefficients.is_empty()) { - op_problem.variable_lower_bounds.resize(op_problem.objective_coefficients.size(), - op_problem.handle_ptr->get_stream()); - thrust::fill(op_problem.handle_ptr->get_thrust_policy(), - op_problem.variable_lower_bounds.begin(), - op_problem.variable_lower_bounds.end(), - f_t(0)); - } - if (op_problem.variable_upper_bounds.is_empty() && - !op_problem.objective_coefficients.is_empty()) { - op_problem.variable_upper_bounds.resize(op_problem.objective_coefficients.size(), - op_problem.handle_ptr->get_stream()); - thrust::fill(op_problem.handle_ptr->get_thrust_policy(), - op_problem.variable_upper_bounds.begin(), - op_problem.variable_upper_bounds.end(), - std::numeric_limits::infinity()); - } + set_variable_bounds(op_problem); if (op_problem.variable_types.is_empty() && !op_problem.objective_coefficients.is_empty()) { op_problem.variable_types.resize(op_problem.objective_coefficients.size(), op_problem.handle_ptr->get_stream()); @@ -261,11 +273,11 @@ static bool check_var_bounds_sanity(const detail::problem_t& problem) bool crossing_bounds_detected = thrust::any_of(problem.handle_ptr->get_thrust_policy(), thrust::counting_iterator(0), - thrust::counting_iterator((i_t)problem.variable_lower_bounds.size()), + thrust::counting_iterator((i_t)problem.variable_bounds.size()), [tolerance = problem.tolerances.presolve_absolute_tolerance, - lb = make_span(problem.variable_lower_bounds), - ub = make_span(problem.variable_upper_bounds)] __device__(i_t index) { - return (lb[index] > ub[index] + tolerance); + var_bnd = make_span(problem.variable_bounds)] __device__(i_t index) { + auto var_bounds = var_bnd[index]; + return (get_lower(var_bounds) > get_upper(var_bounds) + tolerance); }); return !crossing_bounds_detected; } @@ -292,12 +304,12 @@ static void round_bounds(detail::problem_t& problem) thrust::for_each(problem.handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), thrust::make_counting_iterator(problem.n_variables), - [lb = make_span(problem.variable_lower_bounds), - ub = make_span(problem.variable_upper_bounds), - types = make_span(problem.variable_types)] __device__(i_t index) { + [bounds = make_span(problem.variable_bounds), + types = make_span(problem.variable_types)] __device__(i_t index) { if (types[index] == var_t::INTEGER) { - lb[index] = ceil(lb[index]); - ub[index] = floor(ub[index]); + using f_t2 = typename type_2::type; + auto bnd = bounds[index]; + bounds[index] = f_t2{ceil(get_lower(bnd)), floor(get_upper(bnd))}; } }); } diff --git a/cpp/src/mip/problem/write_mps.cu b/cpp/src/mip/problem/write_mps.cu index cca3cd5b1e..ce3e562383 100644 --- a/cpp/src/mip/problem/write_mps.cu +++ b/cpp/src/mip/problem/write_mps.cu @@ -36,8 +36,7 @@ void problem_t::write_as_mps(const std::string& path) auto h_reverse_constraints = cuopt::host_copy(reverse_constraints, handle_ptr->get_stream()); auto h_reverse_offsets = cuopt::host_copy(reverse_offsets, handle_ptr->get_stream()); auto h_obj_coeffs = cuopt::host_copy(objective_coefficients, handle_ptr->get_stream()); - auto h_var_lb = cuopt::host_copy(variable_lower_bounds, handle_ptr->get_stream()); - auto h_var_ub = cuopt::host_copy(variable_upper_bounds, handle_ptr->get_stream()); + auto [h_var_lb, h_var_ub] = extract_host_bounds(variable_bounds, handle_ptr); auto h_cstr_lb = cuopt::host_copy(constraint_lower_bounds, handle_ptr->get_stream()); auto h_cstr_ub = cuopt::host_copy(constraint_upper_bounds, handle_ptr->get_stream()); auto h_var_types = cuopt::host_copy(variable_types, handle_ptr->get_stream()); diff --git a/cpp/src/mip/solution/feasibility_test.cuh b/cpp/src/mip/solution/feasibility_test.cuh index eafb8a67cb..5c452e45fd 100644 --- a/cpp/src/mip/solution/feasibility_test.cuh +++ b/cpp/src/mip/solution/feasibility_test.cuh @@ -55,8 +55,8 @@ __global__ void test_variable_bounds_kernel(typename solution_t::view_ printf("inf var %d val %f l %f u %f integer %d\n", v, val, - sol.problem.variable_lower_bounds[v], - sol.problem.variable_upper_bounds[v], + get_lower(sol.problem.variable_bounds[v]), + get_upper(sol.problem.variable_bounds[v]), sol.problem.is_integer_var(v)); } cuopt_assert(isfinite(val), "assignment should be finite!"); @@ -72,8 +72,8 @@ __global__ void test_variable_bounds_kernel(typename solution_t::view_ printf("oob var %d val %f l %f u %f integer %d\n", v, val, - sol.problem.variable_lower_bounds[v], - sol.problem.variable_upper_bounds[v], + get_lower(sol.problem.variable_bounds[v]), + get_upper(sol.problem.variable_bounds[v]), sol.problem.is_integer_var(v)); } cuopt_assert(feasible, "Variables should be feasible"); diff --git a/cpp/src/mip/solution/solution.cu b/cpp/src/mip/solution/solution.cu index cabb0edda0..7c95aaeed2 100644 --- a/cpp/src/mip/solution/solution.cu +++ b/cpp/src/mip/solution/solution.cu @@ -35,11 +35,25 @@ namespace cuopt::linear_programming::detail { +template +rmm::device_uvector get_lower_bounds( + rmm::device_uvector::type> const& bounds, const raft::handle_t* handle_ptr) +{ + using f_t2 = typename type_2::type; + rmm::device_uvector lower_bounds(bounds.size(), handle_ptr->get_stream()); + thrust::transform(handle_ptr->get_thrust_policy(), + bounds.begin(), + bounds.end(), + lower_bounds.begin(), + [] __device__(auto bnd) { return bnd.x; }); + return lower_bounds; +} + template solution_t::solution_t(problem_t& problem_) : problem_ptr(&problem_), handle_ptr(problem_.handle_ptr), - assignment(problem_.variable_lower_bounds, handle_ptr->get_stream()), + assignment(std::move(get_lower_bounds(problem_.variable_bounds, handle_ptr))), lower_excess(problem_.n_constraints, handle_ptr->get_stream()), upper_excess(problem_.n_constraints, handle_ptr->get_stream()), lower_slack(problem_.n_constraints, handle_ptr->get_stream()), @@ -220,16 +234,16 @@ void solution_t::assign_random_within_bounds(f_t ratio_of_vars_to_rand std::vector h_assignment = host_copy(assignment); std::uniform_real_distribution unif_prob(0, 1); - auto variable_lower_bounds = cuopt::host_copy(problem_ptr->variable_lower_bounds); - auto variable_upper_bounds = cuopt::host_copy(problem_ptr->variable_upper_bounds); - auto variable_types = cuopt::host_copy(problem_ptr->variable_types); + auto variable_bounds = cuopt::host_copy(problem_ptr->variable_bounds); + auto variable_types = cuopt::host_copy(problem_ptr->variable_types); problem_ptr->handle_ptr->sync_stream(); - for (size_t i = 0; i < problem_ptr->variable_lower_bounds.size(); ++i) { + for (size_t i = 0; i < problem_ptr->variable_bounds.size(); ++i) { if (only_integers && variable_types[i] != var_t::INTEGER) { continue; } bool skip = unif_prob(rng) > ratio_of_vars_to_random_assign; if (skip) { continue; } - f_t lower_bound = variable_lower_bounds[i]; - f_t upper_bound = variable_upper_bounds[i]; + auto var_bounds = variable_bounds[i]; + auto lower_bound = get_lower(var_bounds); + auto upper_bound = get_upper(var_bounds); if (lower_bound == -std::numeric_limits::infinity()) { h_assignment[i] = upper_bound; } else if (upper_bound == std::numeric_limits::infinity()) { @@ -322,7 +336,7 @@ template void solution_t::compute_objective() { h_obj = compute_objective_from_vec( - assignment, problem_ptr->objective_coefficients, handle_ptr->get_stream()); + assignment, problem_ptr->objective_coefficients, handle_ptr); // to save from memory transactions, don't update the device objective // when needed we can update the device objective here h_user_obj = problem_ptr->get_user_obj_from_solver_obj(h_obj); @@ -443,10 +457,11 @@ i_t solution_t::calculate_similarity_radius(solution_t& othe problem_ptr->integer_indices.end(), cuda::proclaim_return_type( [other_ptr, curr_assignment, p_view = problem_ptr->view()] __device__(i_t idx) -> bool { + auto var_bounds = p_view.variable_bounds[idx]; return diverse_equal(other_ptr[idx], curr_assignment[idx], - p_view.variable_lower_bounds[idx], - p_view.variable_upper_bounds[idx], + get_lower(var_bounds), + get_upper(var_bounds), p_view.is_integer_var(idx), p_view.tolerances.integrality_tolerance); })); @@ -542,17 +557,15 @@ template f_t solution_t::compute_max_variable_violation() { cuopt_assert(problem_ptr->n_variables == assignment.size(), "Size mismatch"); - cuopt_assert(problem_ptr->n_variables == problem_ptr->variable_lower_bounds.size(), - "Size mismatch"); - cuopt_assert(problem_ptr->n_variables == problem_ptr->variable_upper_bounds.size(), - "Size mismatch"); + cuopt_assert(problem_ptr->n_variables == problem_ptr->variable_bounds.size(), "Size mismatch"); return thrust::transform_reduce( handle_ptr->get_thrust_policy(), thrust::make_counting_iterator(0), thrust::make_counting_iterator(0) + problem_ptr->n_variables, cuda::proclaim_return_type([v = view()] __device__(i_t idx) -> f_t { - f_t lower_vio = max(0., v.problem.variable_lower_bounds[idx] - v.assignment[idx]); - f_t upper_vio = max(0., v.assignment[idx] - v.problem.variable_upper_bounds[idx]); + auto var_bounds = v.problem.variable_bounds[idx]; + f_t lower_vio = max(0., get_lower(var_bounds) - v.assignment[idx]); + f_t upper_vio = max(0., v.assignment[idx] - get_upper(var_bounds)); return max(lower_vio, upper_vio); }), 0., diff --git a/cpp/src/mip/solve.cu b/cpp/src/mip/solve.cu index d16de23f28..486544eba3 100644 --- a/cpp/src/mip/solve.cu +++ b/cpp/src/mip/solve.cu @@ -82,9 +82,9 @@ mip_solution_t run_mip(detail::problem_t& problem, thrust::make_counting_iterator(0), thrust::make_counting_iterator(problem.n_variables), [sol = solution.assignment.data(), pb = problem.view()] __device__(i_t index) { - sol[index] = pb.objective_coefficients[index] > 0 - ? pb.variable_lower_bounds[index] - : pb.variable_upper_bounds[index]; + auto bounds = pb.variable_bounds[index]; + sol[index] = pb.objective_coefficients[index] > 0 ? get_lower(bounds) + : get_upper(bounds); }); problem.post_process_solution(solution); solution.compute_objective(); // just to ensure h_user_obj is set diff --git a/cpp/src/mip/utils.cuh b/cpp/src/mip/utils.cuh index 18759737d5..47f1bbc48b 100644 --- a/cpp/src/mip/utils.cuh +++ b/cpp/src/mip/utils.cuh @@ -218,6 +218,20 @@ bool check_integer_equal_on_indices(const rmm::device_uvector& indices, }); } +template +f_t compute_objective_from_vec(const rmm::device_uvector& assignment, + const rmm::device_uvector& objective_coefficients, + const raft::handle_t* handle_ptr) +{ + cuopt_assert(assignment.size() == objective_coefficients.size(), "Size mismatch!"); + f_t computed_obj = thrust::inner_product(handle_ptr->get_thrust_policy(), + assignment.begin(), + assignment.end(), + objective_coefficients.begin(), + 0.); + return computed_obj; +} + template f_t compute_objective_from_vec(const rmm::device_uvector& assignment, const rmm::device_uvector& objective_coefficients, @@ -239,18 +253,18 @@ void clamp_within_var_bounds(rmm::device_uvector& assignment, { cuopt_assert(assignment.size() == problem_ptr->n_variables, "Size mismatch!"); f_t* assignment_ptr = assignment.data(); - thrust::for_each(handle_ptr->get_thrust_policy(), - thrust::make_counting_iterator(0), - thrust::make_counting_iterator(0) + problem_ptr->n_variables, - [assignment_ptr, - lower_bound = problem_ptr->variable_lower_bounds.data(), - upper_bound = problem_ptr->variable_upper_bounds.data()] __device__(i_t idx) { - if (assignment_ptr[idx] < lower_bound[idx]) { - assignment_ptr[idx] = lower_bound[idx]; - } else if (assignment_ptr[idx] > upper_bound[idx]) { - assignment_ptr[idx] = upper_bound[idx]; - } - }); + thrust::for_each( + handle_ptr->get_thrust_policy(), + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(0) + problem_ptr->n_variables, + [assignment_ptr, variable_bound = problem_ptr->variable_bounds.data()] __device__(i_t idx) { + auto bound = variable_bound[idx]; + if (assignment_ptr[idx] < get_lower(bound)) { + assignment_ptr[idx] = get_lower(bound); + } else if (assignment_ptr[idx] > get_upper(bound)) { + assignment_ptr[idx] = get_upper(bound); + } + }); } template diff --git a/cpp/src/utilities/copy_helpers.hpp b/cpp/src/utilities/copy_helpers.hpp index 78593d3e1a..7207ad72aa 100644 --- a/cpp/src/utilities/copy_helpers.hpp +++ b/cpp/src/utilities/copy_helpers.hpp @@ -18,6 +18,7 @@ #pragma once #include +#include #include #include @@ -27,6 +28,93 @@ #include namespace cuopt { + +template +struct type_2 { + using type = void; +}; + +template <> +struct type_2 { + using type = int2; +}; + +template <> +struct type_2 { + using type = float2; +}; + +template <> +struct type_2 { + using type = double2; +}; + +template +struct scalar_type { + using type = void; +}; + +template <> +struct scalar_type { + using type = int; +}; + +template <> +struct scalar_type { + using type = float; +}; + +template <> +struct scalar_type { + using type = double; +}; + +template <> +struct scalar_type { + using type = const int; +}; + +template <> +struct scalar_type { + using type = const float; +}; + +template <> +struct scalar_type { + using type = const double; +}; + +template +raft::device_span::type> make_span_2(rmm::device_uvector& container) +{ + using T2 = typename type_2::type; + static_assert(sizeof(T2) == 2 * sizeof(T)); + return raft::device_span(reinterpret_cast(container.data()), + sizeof(T) * container.size() / sizeof(T2)); +} + +template +raft::device_span::type> make_span_2( + rmm::device_uvector const& container) +{ + using T2 = typename type_2::type; + static_assert(sizeof(T2) == 2 * sizeof(T)); + return raft::device_span(reinterpret_cast(container.data()), + sizeof(T) * container.size() / sizeof(T2)); +} + +template +__host__ __device__ inline typename scalar_type::type& get_lower(f_t2& val) +{ + return val.x; +} + +template +__host__ __device__ inline typename scalar_type::type& get_upper(f_t2& val) +{ + return val.y; +} + /** * @brief Simple utility function to copy device ptr to host * @@ -253,4 +341,22 @@ inline void expand_device_copy(rmm::device_uvector& dst_vec, raft::copy(dst_vec.data(), src_vec.data(), src_vec.size(), stream_view); } +template +std::tuple, std::vector> extract_host_bounds( + const rmm::device_uvector& variable_bounds, const raft::handle_t* handle_ptr) +{ + rmm::device_uvector var_lb(variable_bounds.size(), handle_ptr->get_stream()); + rmm::device_uvector var_ub(variable_bounds.size(), handle_ptr->get_stream()); + thrust::transform( + handle_ptr->get_thrust_policy(), + variable_bounds.begin(), + variable_bounds.end(), + thrust::make_zip_iterator(thrust::make_tuple(var_lb.begin(), var_ub.begin())), + [] __device__(auto i) { return thrust::make_tuple(get_lower(i), get_upper(i)); }); + handle_ptr->sync_stream(); + auto h_var_lb = cuopt::host_copy(var_lb); + auto h_var_ub = cuopt::host_copy(var_ub); + return std::make_tuple(h_var_lb, h_var_ub); +} + } // namespace cuopt diff --git a/cpp/tests/mip/CMakeLists.txt b/cpp/tests/mip/CMakeLists.txt index 020c537f6a..b9fd249a56 100644 --- a/cpp/tests/mip/CMakeLists.txt +++ b/cpp/tests/mip/CMakeLists.txt @@ -27,9 +27,6 @@ ConfigureTest(ELIM_VAR_REMAP_TEST ConfigureTest(STANDARDIZATION_TEST ${CMAKE_CURRENT_SOURCE_DIR}/bounds_standardization_test.cu ) -ConfigureTest(LB_TEST - ${CMAKE_CURRENT_SOURCE_DIR}/load_balancing_test.cu -) ConfigureTest(MULTI_PROBE_TEST ${CMAKE_CURRENT_SOURCE_DIR}/multi_probe_test.cu ) diff --git a/cpp/tests/mip/elim_var_remap_test.cu b/cpp/tests/mip/elim_var_remap_test.cu index aeb48fe5d6..c486d98c81 100644 --- a/cpp/tests/mip/elim_var_remap_test.cu +++ b/cpp/tests/mip/elim_var_remap_test.cu @@ -100,8 +100,8 @@ void test_elim_var_remap(std::string test_instance) auto fixed_vars = select_k_random(problem.n_variables - 1, 5); for (auto& v : fixed_vars) { double v_val = -v - 1; - problem.variable_lower_bounds.set_element(v, v_val, handle_.get_stream()); - problem.variable_upper_bounds.set_element(v, v_val, handle_.get_stream()); + double2 val = double2{v_val, v_val}; + problem.variable_bounds.set_element(v, val, handle_.get_stream()); full_assignment.set_element(v, v_val, handle_.get_stream()); } // Set free var assignments to 0 @@ -182,8 +182,8 @@ void test_elim_var_solution(std::string test_instance) auto fixed_vars = select_k_random(standardized_problem.n_variables - 1, 5); for (auto& v : fixed_vars) { double v_val = opt_sol_1.get_solution().element(v, handle_.get_stream()); - sub_problem.variable_lower_bounds.set_element(v, v_val, handle_.get_stream()); - sub_problem.variable_upper_bounds.set_element(v, v_val, handle_.get_stream()); + double2 val = double2{v_val, v_val}; + sub_problem.variable_bounds.set_element(v, val, handle_.get_stream()); } handle_.sync_stream(); diff --git a/cpp/tests/mip/load_balancing_test.cu b/cpp/tests/mip/load_balancing_test.cu index deed9ea85a..fb0d8b6e86 100644 --- a/cpp/tests/mip/load_balancing_test.cu +++ b/cpp/tests/mip/load_balancing_test.cu @@ -59,15 +59,14 @@ std::tuple, std::vector, std::vector> select_k_ auto seed = std::random_device{}(); std::cerr << "Tested with seed " << seed << "\n"; problem.compute_n_integer_vars(); - auto v_lb = host_copy(problem.variable_lower_bounds); - auto v_ub = host_copy(problem.variable_upper_bounds); + auto v_bnd = host_copy(problem.variable_bounds); auto int_var_id = host_copy(problem.integer_indices); - int_var_id.erase(std::remove_if(int_var_id.begin(), - int_var_id.end(), - [v_lb, v_ub](auto id) { - return !(std::isfinite(v_lb[id]) && std::isfinite(v_ub[id])); - }), - int_var_id.end()); + int_var_id.erase( + std::remove_if( + int_var_id.begin(), + int_var_id.end(), + [v_bnd](auto id) { return !(std::isfinite(v_bnd[id].x) && std::isfinite(v_bnd[id].y)); }), + int_var_id.end()); sample_size = std::min(sample_size, static_cast(int_var_id.size())); std::vector random_int_vars; std::mt19937 m{seed}; @@ -77,11 +76,11 @@ std::tuple, std::vector, std::vector> select_k_ std::vector probe_1(sample_size); for (int i = 0; i < static_cast(random_int_vars.size()); ++i) { if (i % 2) { - probe_0[i] = v_lb[random_int_vars[i]]; - probe_1[i] = v_ub[random_int_vars[i]]; + probe_0[i] = v_bnd[random_int_vars[i]].x; + probe_1[i] = v_bnd[random_int_vars[i]].y; } else { - probe_1[i] = v_lb[random_int_vars[i]]; - probe_0[i] = v_ub[random_int_vars[i]]; + probe_1[i] = v_bnd[random_int_vars[i]].x; + probe_0[i] = v_bnd[random_int_vars[i]].y; } } return std::make_tuple(std::move(random_int_vars), std::move(probe_0), std::move(probe_1)); diff --git a/cpp/tests/mip/multi_probe_test.cu b/cpp/tests/mip/multi_probe_test.cu index 63cf93c792..1473c84bff 100644 --- a/cpp/tests/mip/multi_probe_test.cu +++ b/cpp/tests/mip/multi_probe_test.cu @@ -58,15 +58,15 @@ std::tuple, std::vector, std::vector> select_k_ auto seed = std::random_device{}(); std::cerr << "Tested with seed " << seed << "\n"; problem.compute_n_integer_vars(); - auto v_lb = host_copy(problem.variable_lower_bounds); - auto v_ub = host_copy(problem.variable_upper_bounds); - auto int_var_id = host_copy(problem.integer_indices); - int_var_id.erase(std::remove_if(int_var_id.begin(), - int_var_id.end(), - [v_lb, v_ub](auto id) { - return !(std::isfinite(v_lb[id]) && std::isfinite(v_ub[id])); - }), - int_var_id.end()); + auto [v_lb, v_ub] = extract_host_bounds(problem.variable_bounds, problem.handle_ptr); + auto int_var_id = host_copy(problem.integer_indices); + int_var_id.erase( + std::remove_if(int_var_id.begin(), + int_var_id.end(), + [v_lb_sp = v_lb, v_ub_sp = v_ub](auto id) { + return !(std::isfinite(v_lb_sp[id]) && std::isfinite(v_ub_sp[id])); + }), + int_var_id.end()); sample_size = std::min(sample_size, static_cast(int_var_id.size())); std::vector random_int_vars; std::mt19937 m{seed};