diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index f6b30e72c..64b1bb4bb 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -17,6 +17,7 @@ #include "initial_solution_reader.hpp" #include "mip_test_instances.hpp" +#include #include #include #include @@ -24,7 +25,6 @@ #include #include -#include #include #include diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 2986de018..3d4192105 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -15,6 +15,7 @@ * limitations under the License. */ +#include #include #include @@ -30,8 +31,10 @@ #include #include +#include +#include #include -#include +#include namespace cuopt::linear_programming::dual_simplex { @@ -123,7 +126,7 @@ bool check_guess(const lp_problem_t& original_lp, } template -void set_uninitialized_steepest_edge_norms(i_t n, std::vector& edge_norms) +void set_uninitialized_steepest_edge_norms(std::vector& edge_norms) { for (i_t j = 0; j < edge_norms.size(); ++j) { if (edge_norms[j] <= 0.0) { edge_norms[j] = 1e-4; } @@ -184,17 +187,6 @@ dual::status_t convert_lp_status_to_dual_status(lp_status_t status) } } -} // namespace - -template -f_t branch_and_bound_t::get_upper_bound() -{ - mutex_upper.lock(); - const f_t upper_bound = upper_bound_; - mutex_upper.unlock(); - return upper_bound; -} - template f_t sgn(f_t x) { @@ -207,7 +199,8 @@ f_t relative_gap(f_t obj_value, f_t lower_bound) f_t user_mip_gap = obj_value == 0.0 ? (lower_bound == 0.0 ? 0.0 : std::numeric_limits::infinity()) : std::abs(obj_value - lower_bound) / std::abs(obj_value); - if (user_mip_gap != user_mip_gap) { return std::numeric_limits::infinity(); } + // Handle NaNs (i.e., NaN != NaN) + if (std::isnan(user_mip_gap)) { return std::numeric_limits::infinity(); } return user_mip_gap; } @@ -225,33 +218,81 @@ std::string user_mip_gap(f_t obj_value, f_t lower_bound) } } +} // namespace + +template +branch_and_bound_t::branch_and_bound_t( + const user_problem_t& user_problem, + const simplex_solver_settings_t& solver_settings) + : original_problem_(user_problem), + settings_(solver_settings), + original_lp_(1, 1, 1), + incumbent_(1), + root_relax_soln_(1, 1), + pc_(1) +{ + stats_.start_time = tic(); + convert_user_problem(original_problem_, settings_, original_lp_, new_slacks_); + full_variable_types(original_problem_, original_lp_, var_types_); + + mutex_upper_.lock(); + upper_bound_ = inf; + mutex_upper_.unlock(); + + mutex_lower_.lock(); + lower_bound_ = -inf; + mutex_lower_.unlock(); + + mutex_branching_.lock(); + currently_branching_ = false; + mutex_branching_.unlock(); +} + +template +f_t branch_and_bound_t::get_upper_bound() +{ + mutex_upper_.lock(); + const f_t upper_bound = upper_bound_; + mutex_upper_.unlock(); + return upper_bound; +} + +template +f_t branch_and_bound_t::get_lower_bound() +{ + mutex_lower_.lock(); + const f_t lower_bound = lower_bound_; + mutex_lower_.unlock(); + return lower_bound; +} + template void branch_and_bound_t::set_new_solution(const std::vector& solution) { - if (solution.size() != original_problem.num_cols) { - settings.log.printf( - "Solution size mismatch %ld %d\n", solution.size(), original_problem.num_cols); + if (solution.size() != original_problem_.num_cols) { + settings_.log.printf( + "Solution size mismatch %ld %d\n", solution.size(), original_problem_.num_cols); } std::vector crushed_solution; crush_primal_solution( - original_problem, original_lp, solution, new_slacks, crushed_solution); - f_t obj = compute_objective(original_lp, crushed_solution); + original_problem_, original_lp_, solution, new_slacks_, crushed_solution); + f_t obj = compute_objective(original_lp_, crushed_solution); bool is_feasible = false; bool attempt_repair = false; - mutex_upper.lock(); + mutex_upper_.lock(); if (obj < upper_bound_) { f_t primal_err; f_t bound_err; i_t num_fractional; is_feasible = check_guess( - original_lp, settings, var_types, crushed_solution, primal_err, bound_err, num_fractional); + original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); if (is_feasible) { upper_bound_ = obj; } else { attempt_repair = true; constexpr bool verbose = false; if (verbose) { - settings.log.printf( + settings_.log.printf( "Injected solution infeasible. Constraint error %e bound error %e integer infeasible " "%d\n", primal_err, @@ -260,68 +301,65 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu } } } - mutex_upper.unlock(); + mutex_upper_.unlock(); if (is_feasible) { - mutex_lower.lock(); - f_t lower_bound = lower_bound_; - mutex_lower.unlock(); - mutex_branching.lock(); - bool currently_branching = currently_branching; - mutex_branching.unlock(); + mutex_branching_.lock(); + bool currently_branching = currently_branching_; + mutex_branching_.unlock(); if (currently_branching) { - settings.log.printf( + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap = user_mip_gap(user_obj, user_lower); + + settings_.log.printf( "H %+13.6e %+10.6e %s %9.2f\n", - compute_user_objective(original_lp, obj), - compute_user_objective(original_lp, lower_bound), - user_mip_gap(compute_user_objective(original_lp, obj), - compute_user_objective(original_lp, lower_bound)) - .c_str(), - toc(start_time)); + user_obj, + user_lower, + gap.c_str(), + toc(stats_.start_time)); } else { - settings.log.printf("New solution from primal heuristics. Objective %+.6e. Time %.2f\n", - compute_user_objective(original_lp, obj), - toc(start_time)); + settings_.log.printf("New solution from primal heuristics. Objective %+.6e. Time %.2f\n", + compute_user_objective(original_lp_, obj), + toc(stats_.start_time)); } } if (attempt_repair) { - mutex_repair.lock(); - repair_queue.push_back(crushed_solution); - mutex_repair.unlock(); + mutex_repair_.lock(); + repair_queue_.push_back(crushed_solution); + mutex_repair_.unlock(); } } template -bool branch_and_bound_t::repair_solution( - const std::vector& root_vstatus, - const std::vector& edge_norms, - const std::vector& potential_solution, - f_t& repaired_obj, - std::vector& repaired_solution) const +bool branch_and_bound_t::repair_solution(const std::vector& edge_norms, + const std::vector& potential_solution, + f_t& repaired_obj, + std::vector& repaired_solution) const { bool feasible = false; repaired_obj = std::numeric_limits::quiet_NaN(); - i_t n = original_lp.num_cols; + i_t n = original_lp_.num_cols; assert(potential_solution.size() == n); - lp_problem_t repair_lp = original_lp; + lp_problem_t repair_lp = original_lp_; // Fix integer variables for (i_t j = 0; j < n; ++j) { - if (var_types[j] == variable_type_t::INTEGER) { + if (var_types_[j] == variable_type_t::INTEGER) { const f_t fixed_val = std::round(potential_solution[j]); repair_lp.lower[j] = fixed_val; repair_lp.upper[j] = fixed_val; } } - lp_solution_t lp_solution(original_lp.num_rows, original_lp.num_cols); + lp_solution_t lp_solution(original_lp_.num_rows, original_lp_.num_cols); i_t iter = 0; f_t lp_start_time = tic(); - simplex_solver_settings_t lp_settings = settings; - std::vector vstatus = root_vstatus; + simplex_solver_settings_t lp_settings = settings_; + std::vector vstatus = root_vstatus_; lp_settings.set_log(false); lp_settings.inside_mip = true; std::vector leaf_edge_norms = edge_norms; @@ -334,12 +372,17 @@ bool branch_and_bound_t::repair_solution( f_t primal_error; f_t bound_error; i_t num_fractional; - feasible = check_guess( - original_lp, settings, var_types, lp_solution.x, primal_error, bound_error, num_fractional); - repaired_obj = compute_objective(original_lp, repaired_solution); + feasible = check_guess(original_lp_, + settings_, + var_types_, + lp_solution.x, + primal_error, + bound_error, + num_fractional); + repaired_obj = compute_objective(original_lp_, repaired_solution); constexpr bool verbose = false; if (verbose) { - settings.log.printf( + settings_.log.printf( "After repair: feasible %d primal error %e bound error %e fractional %d. Objective %e\n", feasible, primal_error, @@ -353,62 +396,501 @@ bool branch_and_bound_t::repair_solution( } template -branch_and_bound_t::branch_and_bound_t( - const user_problem_t& user_problem, - const simplex_solver_settings_t& solver_settings) - : original_problem(user_problem), settings(solver_settings), original_lp(1, 1, 1), incumbent(1) +void branch_and_bound_t::repair_heuristic_solutions() { - start_time = tic(); - convert_user_problem(original_problem, settings, original_lp, new_slacks); - full_variable_types(original_problem, original_lp, var_types); + // Check if there are any solutions to repair + std::vector> to_repair; + mutex_repair_.lock(); + if (repair_queue_.size() > 0) { + to_repair = repair_queue_; + repair_queue_.clear(); + } + mutex_repair_.unlock(); + + if (to_repair.size() > 0) { + settings_.log.debug("Attempting to repair %ld injected solutions\n", to_repair.size()); + for (const std::vector& potential_solution : to_repair) { + std::vector repaired_solution; + f_t repaired_obj; + bool is_feasible = + repair_solution(edge_norms_, potential_solution, repaired_obj, repaired_solution); + if (is_feasible) { + mutex_upper_.lock(); + + if (repaired_obj < upper_bound_) { + upper_bound_ = repaired_obj; + incumbent_.set_incumbent_solution(repaired_obj, repaired_solution); + + f_t obj = compute_user_objective(original_lp_, repaired_obj); + f_t lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string user_gap = user_mip_gap(obj, lower); + settings_.log.printf( + "H %+13.6e %+10.6e %s %9.2f\n", + obj, + lower, + user_gap.c_str(), + toc(stats_.start_time)); + + if (settings_.solution_callback != nullptr) { + std::vector original_x; + uncrush_primal_solution(original_problem_, original_lp_, repaired_solution, original_x); + settings_.solution_callback(original_x, repaired_obj); + } + } - mutex_upper.lock(); - upper_bound_ = inf; - mutex_upper.unlock(); + mutex_upper_.unlock(); + } + } + } +} - mutex_lower.lock(); - lower_bound_ = -inf; - mutex_lower.unlock(); +template +void branch_and_bound_t::branch(mip_node_t* parent_node, + i_t branch_var, + f_t branch_var_val, + const std::vector& parent_vstatus) +{ + // down child + auto down_child = std::make_unique>( + original_lp_, parent_node, ++stats_.num_nodes, branch_var, 0, branch_var_val, parent_vstatus); + + graphviz_edge( + settings_, parent_node, down_child.get(), branch_var, 0, std::floor(branch_var_val)); + + // up child + auto up_child = std::make_unique>( + original_lp_, parent_node, ++stats_.num_nodes, branch_var, 1, branch_var_val, parent_vstatus); + + graphviz_edge(settings_, parent_node, up_child.get(), branch_var, 1, std::ceil(branch_var_val)); + + assert(parent_vstatus.size() == original_lp_.num_cols); + parent_node->add_children(std::move(down_child), + std::move(up_child)); // child pointers moved into the tree +} + +template +void branch_and_bound_t::update_tree(mip_node_t* node_ptr, node_status_t status) +{ + std::vector*> stack; + node_ptr->set_status(status, stack); + remove_fathomed_nodes(stack); +} + +template +void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, + const std::vector& leaf_solution, + i_t leaf_depth, + char symbol) +{ + bool send_solution = false; + i_t nodes_explored = stats_.nodes_explored; + i_t nodes_unexplored = stats_.nodes_unexplored; + f_t gap; + + mutex_upper_.lock(); + if (leaf_objective < upper_bound_) { + incumbent_.set_incumbent_solution(leaf_objective, leaf_solution); + upper_bound_ = leaf_objective; + f_t lower_bound = get_lower_bound(); + gap = upper_bound_ - lower_bound; + f_t obj = compute_user_objective(original_lp_, upper_bound_); + f_t lower = compute_user_objective(original_lp_, lower_bound); + settings_.log.printf("%c%8d %8lu %+13.6e %+10.6e %4d %7.1e %s %9.2f\n", + symbol, + nodes_explored, + nodes_unexplored, + obj, + lower, + leaf_depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + user_mip_gap(obj, lower).c_str(), + toc(stats_.start_time)); + + send_solution = true; + } + + if (send_solution && settings_.solution_callback != nullptr) { + std::vector original_x; + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, original_x); + settings_.solution_callback(original_x, upper_bound_); + } + mutex_upper_.unlock(); + + if (send_solution) { + mutex_gap_.lock(); + gap_ = gap; + mutex_gap_.unlock(); + } +} + +template +dual::status_t branch_and_bound_t::node_dual_simplex( + i_t leaf_id, + lp_problem_t& leaf_problem, + std::vector& leaf_vstatus, + lp_solution_t& leaf_solution, + std::vector& bounds_changed, + csc_matrix_t& Arow, + f_t upper_bound, + logger_t& log) +{ + i_t node_iter = 0; + assert(leaf_vstatus.size() == leaf_problem.num_cols); + f_t lp_start_time = tic(); + std::vector leaf_edge_norms = edge_norms_; // = node.steepest_edge_norms; + simplex_solver_settings_t lp_settings = settings_; + lp_settings.set_log(false); + lp_settings.cut_off = upper_bound + settings_.dual_tol; + lp_settings.inside_mip = 2; + + // in B&B we only have equality constraints, leave it empty for default + std::vector row_sense; + bool feasible = + bound_strengthening(row_sense, lp_settings, leaf_problem, Arow, var_types_, bounds_changed); + + dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; + + if (feasible) { + lp_status = dual_phase2(2, + 0, + lp_start_time, + leaf_problem, + lp_settings, + leaf_vstatus, + leaf_solution, + node_iter, + leaf_edge_norms); + + if (lp_status == dual::status_t::NUMERICAL) { + log.printf("Numerical issue node %d. Resolving from scratch.\n", leaf_id); + lp_status_t second_status = solve_linear_program_advanced( + leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); + lp_status = convert_lp_status_to_dual_status(second_status); + } + } else { + log.printf("Infeasible after bounds strengthening. Fathoming node %d.\n", leaf_id); + } + + mutex_stats_.lock(); + stats_.total_lp_solve_time += toc(lp_start_time); + stats_.total_lp_iters += node_iter; + mutex_stats_.unlock(); + + return lp_status; +} + +template +mip_status_t branch_and_bound_t::solve_node_lp( + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + const std::vector& var_types, + f_t upper_bound) +{ + logger_t log; + log.log = false; + std::vector& leaf_vstatus = node_ptr->vstatus; + lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); + + // Set the correct bounds for the leaf problem + leaf_problem.lower = original_lp_.lower; + leaf_problem.upper = original_lp_.upper; + + std::vector bounds_changed(original_lp_.num_cols, false); + // Technically, we can get the already strengthened bounds from the node/parent instead of + // getting it from the original problem and re-strengthening. But this requires storing + // two vectors at each node and potentially cause memory issues + node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); + + dual::status_t lp_status = node_dual_simplex(node_ptr->node_id, + leaf_problem, + leaf_vstatus, + leaf_solution, + bounds_changed, + Arow, + upper_bound, + settings_.log); + + if (lp_status == dual::status_t::DUAL_UNBOUNDED) { + node_ptr->lower_bound = inf; + graphviz_node(settings_, node_ptr, "infeasible", 0.0); + update_tree(node_ptr, node_status_t::INFEASIBLE); + // Node was infeasible. Do not branch + + } else if (lp_status == dual::status_t::CUTOFF) { + node_ptr->lower_bound = upper_bound; + f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); + graphviz_node(settings_, node_ptr, "cut off", leaf_objective); + update_tree(node_ptr, node_status_t::FATHOMED); + // Node was cut off. Do not branch + } else if (lp_status == dual::status_t::OPTIMAL) { + // LP was feasible + std::vector fractional; + const i_t leaf_fractional = + fractional_variables(settings_, leaf_solution.x, var_types_, fractional); + f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); + graphviz_node(settings_, node_ptr, "lower bound", leaf_objective); + + mutex_pc_.lock(); + pc_.update_pseudo_costs(node_ptr, leaf_objective); + mutex_pc_.unlock(); + + node_ptr->lower_bound = leaf_objective; + + constexpr f_t fathom_tol = 1e-5; + if (leaf_fractional == 0) { + add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, 'B'); + graphviz_node(settings_, node_ptr, "integer feasible", leaf_objective); + update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); + + } else if (leaf_objective <= upper_bound + fathom_tol) { + // Choose fractional variable to branch on + mutex_pc_.lock(); + const i_t branch_var = pc_.variable_selection( + fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, log); + mutex_pc_.unlock(); + + assert(leaf_vstatus.size() == leaf_problem.num_cols); + branch(node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus); + node_ptr->status = node_status_t::HAS_CHILDREN; + + } else { + graphviz_node(settings_, node_ptr, "fathomed", leaf_objective); + update_tree(node_ptr, node_status_t::FATHOMED); + } + } else { + graphviz_node(settings_, node_ptr, "numerical", 0.0); + settings_.log.printf("Encountered LP status %d. This indicates a numerical issue.\n", + lp_status); + return mip_status_t::NUMERICAL; + } + + return mip_status_t::UNSET; +} + +template +mip_status_t branch_and_bound_t::explore_tree(i_t branch_var, + mip_solution_t& solution) +{ + mip_status_t status = mip_status_t::UNSET; + logger_t log; + log.log = false; + auto compare = [](mip_node_t* a, mip_node_t* b) { + return a->lower_bound > + b->lower_bound; // True if a comes before b, elements that come before are output last + }; + + std::priority_queue*, std::vector*>, decltype(compare)> + heap(compare); + + mip_node_t root_node(root_objective_, root_vstatus_); + graphviz_node(settings_, &root_node, "lower bound", root_objective_); + + branch(&root_node, branch_var, root_relax_soln_.x[branch_var], root_vstatus_); + + // the stack does not own the unique_ptr the tree does + heap.push(root_node.get_down_child()); + heap.push(root_node.get_up_child()); + + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + csc_matrix_t Arow(1, 1, 1); + original_lp_.A.transpose(Arow); + + f_t lower_bound = get_lower_bound(); + f_t gap = get_upper_bound() - lower_bound; + f_t last_log = 0; + + while (gap > settings_.absolute_mip_gap_tol && + relative_gap(get_upper_bound(), lower_bound) > settings_.relative_mip_gap_tol && + heap.size() > 0) { + repair_heuristic_solutions(); + + // Get a node off the heap + mip_node_t* node_ptr = heap.top(); + heap.pop(); + stats_.nodes_explored++; + + f_t upper_bound = get_upper_bound(); + if (upper_bound < node_ptr->lower_bound) { + // This node was put on the heap earlier but its lower bound is now greater than the current + // upper bound + update_tree(node_ptr, node_status_t::FATHOMED); + graphviz_node(settings_, node_ptr, "cutoff", node_ptr->lower_bound); + continue; + } + mutex_lower_.lock(); + lower_bound = lower_bound_ = node_ptr->lower_bound; + mutex_lower_.unlock(); + + mutex_gap_.lock(); + gap_ = gap = upper_bound - lower_bound; + mutex_gap_.unlock(); + + i_t nodes_explored = stats_.nodes_explored; + f_t now = toc(stats_.start_time); + f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); + if ((nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || + nodes_explored < 1000) && + (time_since_log >= 1) || + (time_since_log > 60) || now > settings_.time_limit) { + f_t user_obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, lower_bound); + std::string user_gap = user_mip_gap(user_obj, user_lower); + + settings_.log.printf(" %8d %8lu %+13.6e %+10.6e %4d %7.1e %s %9.2f\n", + nodes_explored, + heap.size(), + user_obj, + user_lower, + node_ptr->depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + user_gap.c_str(), + now); + last_log = tic(); + } + + if (toc(stats_.start_time) > settings_.time_limit) { + settings_.log.printf("Hit time limit. Stopping\n"); + status = mip_status_t::TIME_LIMIT; + break; + } + + status = solve_node_lp(node_ptr, leaf_problem, Arow, var_types_, upper_bound); + + if (status == mip_status_t::NUMERICAL) { break; } + + if (node_ptr->status == node_status_t::HAS_CHILDREN) { + // the heap does not own the unique_ptr the tree does + heap.push(node_ptr->get_down_child()); + heap.push(node_ptr->get_up_child()); + } + } + + stats_.nodes_unexplored = heap.size(); + + if (stats_.nodes_unexplored == 0) { + mutex_lower_.lock(); + lower_bound = lower_bound_ = root_node.lower_bound; + mutex_lower_.unlock(); + + mutex_gap_.lock(); + gap_ = gap = get_upper_bound() - lower_bound; + mutex_gap_.unlock(); + } + + return status; +} + +template +mip_status_t branch_and_bound_t::dive(i_t branch_var, mip_solution_t& solution) +{ + mip_status_t status = mip_status_t::UNSET; + + logger_t log; + log.log = false; - mutex_branching.lock(); - currently_branching = false; - mutex_branching.unlock(); + std::vector*> node_stack; + + mip_node_t root_node(root_objective_, root_vstatus_); + graphviz_node(settings_, &root_node, "lower bound", root_objective_); + + branch(&root_node, branch_var, root_relax_soln_.x[branch_var], root_vstatus_); + + // the stack does not own the unique_ptr the tree does + node_stack.push_back(root_node.get_down_child()); + node_stack.push_back(root_node.get_up_child()); + + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + + csc_matrix_t Arow(1, 1, 1); + original_lp_.A.transpose(Arow); + + f_t lower_bound = get_lower_bound(); + f_t gap = get_upper_bound() - lower_bound; + i_t nodes_explored = 0; + + while (node_stack.size() > 0) { + // Get a node off the stack + mip_node_t* node_ptr = node_stack.back(); + node_stack.pop_back(); + nodes_explored++; + + f_t upper_bound = get_upper_bound(); + lower_bound = get_lower_bound(); + gap = upper_bound - lower_bound; + + if (gap < settings_.absolute_mip_gap_tol && + relative_gap(get_upper_bound(), lower_bound) < settings_.relative_mip_gap_tol) { + update_tree(node_ptr, node_status_t::FATHOMED); + continue; + } + + if (toc(stats_.start_time) > settings_.time_limit) { + status = mip_status_t::TIME_LIMIT; + break; + } + + status = solve_node_lp(node_ptr, leaf_problem, Arow, var_types_, upper_bound); + + if (status == mip_status_t::NUMERICAL) { continue; } + + if (node_ptr->status == node_status_t::HAS_CHILDREN) { + // Martin's child selection + const i_t branch_var = node_ptr->get_down_child()->branch_var; + const f_t branch_var_val = node_ptr->get_down_child()->fractional_val; + const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); + const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); + const f_t down_dist = branch_var_val - down_val; + const f_t up_dist = up_val - branch_var_val; + + if (down_dist < up_dist) { + node_stack.push_back(node_ptr->get_up_child()); + node_stack.push_back(node_ptr->get_down_child()); + } else { + node_stack.push_back(node_ptr->get_down_child()); + node_stack.push_back(node_ptr->get_up_child()); + } + } + } + + return status; } template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { mip_status_t status = mip_status_t::UNSET; - mip_solution_t incumbent(original_lp.num_cols); - if (guess.size() != 0) { + + if (guess_.size() != 0) { std::vector crushed_guess; - crush_primal_solution(original_problem, original_lp, guess, new_slacks, crushed_guess); + crush_primal_solution(original_problem_, original_lp_, guess_, new_slacks_, crushed_guess); f_t primal_err; f_t bound_err; i_t num_fractional; const bool feasible = check_guess( - original_lp, settings, var_types, crushed_guess, primal_err, bound_err, num_fractional); + original_lp_, settings_, var_types_, crushed_guess, primal_err, bound_err, num_fractional); if (feasible) { - const f_t computed_obj = compute_objective(original_lp, crushed_guess); - mutex_upper.lock(); - incumbent.set_incumbent_solution(computed_obj, crushed_guess); + const f_t computed_obj = compute_objective(original_lp_, crushed_guess); + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(computed_obj, crushed_guess); upper_bound_ = computed_obj; - mutex_upper.unlock(); + mutex_upper_.unlock(); } } - lp_solution_t root_relax_soln(original_lp.num_rows, original_lp.num_cols); - std::vector root_vstatus; - std::vector edge_norms; - settings.log.printf("Solving LP root relaxation\n"); - simplex_solver_settings_t lp_settings = settings; + root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols); + + settings_.log.printf("Solving LP root relaxation\n"); + simplex_solver_settings_t lp_settings = settings_; lp_settings.inside_mip = 1; lp_status_t root_status = solve_linear_program_advanced( - original_lp, start_time, lp_settings, root_relax_soln, root_vstatus, edge_norms); - f_t total_lp_solve_time = toc(start_time); - assert(root_vstatus.size() == original_lp.num_cols); + original_lp_, stats_.start_time, lp_settings, root_relax_soln_, root_vstatus_, edge_norms_); + stats_.total_lp_solve_time = toc(stats_.start_time); + assert(root_vstatus_.size() == original_lp_.num_cols); if (root_status == lp_status_t::INFEASIBLE) { - settings.log.printf("MIP Infeasible\n"); + settings_.log.printf("MIP Infeasible\n"); // FIXME: rarely dual simplex detects infeasible whereas it is feasible. // to add a small safety net, check if there is a primal solution already. // Uncomment this if the issue with cost266-UUE is resolved @@ -418,444 +900,148 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::INFEASIBLE; } if (root_status == lp_status_t::UNBOUNDED) { - settings.log.printf("MIP Unbounded\n"); - if (settings.heuristic_preemption_callback != nullptr) { - settings.heuristic_preemption_callback(); + settings_.log.printf("MIP Unbounded\n"); + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); } return mip_status_t::UNBOUNDED; } if (root_status == lp_status_t::TIME_LIMIT) { - settings.log.printf("Hit time limit\n"); + settings_.log.printf("Hit time limit\n"); return mip_status_t::TIME_LIMIT; } - set_uninitialized_steepest_edge_norms(original_lp.num_cols, edge_norms); + set_uninitialized_steepest_edge_norms(edge_norms_); - std::vector fractional; - const i_t num_fractional = - fractional_variables(settings, root_relax_soln.x, var_types, fractional); - const f_t root_objective = compute_objective(original_lp, root_relax_soln.x); - if (settings.set_simplex_solution_callback != nullptr) { + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + if (settings_.set_simplex_solution_callback != nullptr) { std::vector original_x; - uncrush_primal_solution(original_problem, original_lp, root_relax_soln.x, original_x); + uncrush_primal_solution(original_problem_, original_lp_, root_relax_soln_.x, original_x); std::vector original_dual; std::vector original_z; - uncrush_dual_solution(original_problem, - original_lp, - root_relax_soln.y, - root_relax_soln.z, + uncrush_dual_solution(original_problem_, + original_lp_, + root_relax_soln_.y, + root_relax_soln_.z, original_dual, original_z); - settings.set_simplex_solution_callback( - original_x, original_dual, compute_user_objective(original_lp, root_objective)); + settings_.set_simplex_solution_callback( + original_x, original_dual, compute_user_objective(original_lp_, root_objective_)); } - mutex_lower.lock(); - f_t lower_bound = lower_bound_ = root_objective; - mutex_lower.unlock(); + mutex_lower_.lock(); + lower_bound_ = root_objective_; + mutex_lower_.unlock(); + + std::vector fractional; + const i_t num_fractional = + fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); if (num_fractional == 0) { - mutex_upper.lock(); - incumbent.set_incumbent_solution(root_objective, root_relax_soln.x); - upper_bound_ = root_objective; - mutex_upper.unlock(); + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); + upper_bound_ = root_objective_; + mutex_upper_.unlock(); // We should be done here - uncrush_primal_solution(original_problem, original_lp, incumbent.x, solution.x); - solution.objective = incumbent.objective; - solution.lower_bound = lower_bound; + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); + solution.objective = incumbent_.objective; + solution.lower_bound = lower_bound_; solution.nodes_explored = 0; - solution.simplex_iterations = root_relax_soln.iterations; - settings.log.printf("Optimal solution found at root node. Objective %.16e. Time %.2f.\n", - compute_user_objective(original_lp, root_objective), - toc(start_time)); - if (settings.solution_callback != nullptr) { - settings.solution_callback(solution.x, solution.objective); + solution.simplex_iterations = root_relax_soln_.iterations; + settings_.log.printf("Optimal solution found at root node. Objective %.16e. Time %.2f.\n", + compute_user_objective(original_lp_, root_objective_), + toc(stats_.start_time)); + if (settings_.solution_callback != nullptr) { + settings_.solution_callback(solution.x, solution.objective); } - if (settings.heuristic_preemption_callback != nullptr) { - settings.heuristic_preemption_callback(); + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); } return mip_status_t::OPTIMAL; } - pseudo_costs_t pc(original_lp.num_cols); - strong_branching(original_lp, - settings, - start_time, - var_types, - root_relax_soln.x, + pc_.resize(original_lp_.num_cols); + strong_branching(original_lp_, + settings_, + stats_.start_time, + var_types_, + root_relax_soln_.x, fractional, - root_objective, - root_vstatus, - edge_norms, - pc); - - auto compare = [](mip_node_t* a, mip_node_t* b) { - return a->lower_bound > - b->lower_bound; // True if a comes before b, elements that come before are output last - }; - std::priority_queue*, std::vector*>, decltype(compare)> - heap(compare); - i_t num_nodes = 0; - mip_node_t root_node(root_objective, root_vstatus); - graphviz_node(settings, &root_node, "lower bound", lower_bound); + root_objective_, + root_vstatus_, + edge_norms_, + pc_); // Choose variable to branch on logger_t log; - log.log = false; - const i_t branch_var = - pc.variable_selection(fractional, root_relax_soln.x, original_lp.lower, original_lp.upper, log); + log.log = false; + i_t branch_var = pc_.variable_selection( + fractional, root_relax_soln_.x, original_lp_.lower, original_lp_.upper, log); - // down child - std::unique_ptr> down_child = - std::make_unique>(original_lp, - &root_node, - ++num_nodes, - branch_var, - 0, - root_relax_soln.x[branch_var], - root_vstatus); - - graphviz_edge(settings, - &root_node, - down_child.get(), - branch_var, - 0, - std::floor(root_relax_soln.x[branch_var])); - - // up child - std::unique_ptr> up_child = - std::make_unique>(original_lp, - &root_node, - ++num_nodes, - branch_var, - 1, - root_relax_soln.x[branch_var], - root_vstatus); + stats_.total_lp_iters = 0; + stats_.nodes_explored = 0; + stats_.nodes_unexplored = 0; + stats_.num_nodes = 1; - graphviz_edge( - settings, &root_node, up_child.get(), branch_var, 0, std::ceil(root_relax_soln.x[branch_var])); - - assert(root_vstatus.size() == original_lp.num_cols); - heap.push(down_child.get()); // the heap does not own the unique_ptr the tree does - heap.push(up_child.get()); // the heap does not own the unqiue_ptr the tree does - root_node.add_children(std::move(down_child), - std::move(up_child)); // child pointers moved into the tree - lp_problem_t leaf_problem = - original_lp; // Make a copy of the original LP. We will modify its bounds at each leaf - csc_matrix_t Arow(1, 1, 1); - original_lp.A.transpose(Arow); - f_t gap = get_upper_bound() - lower_bound; - i_t nodes_explored = 0; - settings.log.printf( + settings_.log.printf( "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap | " " Time \n"); - mutex_branching.lock(); - currently_branching = true; - mutex_branching.unlock(); - - f_t total_lp_iters = 0.0; - f_t last_log = 0; - while (gap > settings.absolute_mip_gap_tol && - relative_gap(get_upper_bound(), lower_bound) > settings.relative_mip_gap_tol && - heap.size() > 0) { - // Check if there are any solutions to repair - std::vector> to_repair; - mutex_repair.lock(); - if (repair_queue.size() > 0) { - to_repair = repair_queue; - repair_queue.clear(); - } - mutex_repair.unlock(); - - if (to_repair.size() > 0) { - settings.log.debug("Attempting to repair %ld injected solutions\n", to_repair.size()); - for (const std::vector& potential_solution : to_repair) { - std::vector repaired_solution; - f_t repaired_obj; - bool is_feasible = repair_solution( - root_vstatus, edge_norms, potential_solution, repaired_obj, repaired_solution); - if (is_feasible) { - mutex_upper.lock(); - if (repaired_obj < upper_bound_) { - upper_bound_ = repaired_obj; - incumbent.set_incumbent_solution(repaired_obj, repaired_solution); - - settings.log.printf( - "H %+13.6e %+10.6e %s %9.2f\n", - compute_user_objective(original_lp, repaired_obj), - compute_user_objective(original_lp, lower_bound), - user_mip_gap(compute_user_objective(original_lp, repaired_obj), - compute_user_objective(original_lp, lower_bound)) - .c_str(), - toc(start_time)); - if (settings.solution_callback != nullptr) { - std::vector original_x; - uncrush_primal_solution(original_problem, original_lp, repaired_solution, original_x); - settings.solution_callback(original_x, repaired_obj); - } - } - mutex_upper.unlock(); - } - } - } - - // Get a node off the heap - mip_node_t* node_ptr = heap.top(); - heap.pop(); // Remove node from the heap - f_t upper_bound = get_upper_bound(); - if (upper_bound < node_ptr->lower_bound) { - // This node was put on the heap earlier but its lower bound is now greater than the current - // upper bound - std::vector*> stack; - node_ptr->set_status(node_status_t::FATHOMED, stack); - graphviz_node(settings, node_ptr, "cutoff", node_ptr->lower_bound); - remove_fathomed_nodes(stack); - continue; - } - mutex_lower.lock(); - lower_bound_ = lower_bound = node_ptr->lower_bound; - mutex_lower.unlock(); - gap = upper_bound - lower_bound; - const i_t leaf_depth = node_ptr->depth; - f_t now = toc(start_time); - f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); - if ((nodes_explored % 1000 == 0 || gap < 10 * settings.absolute_mip_gap_tol || - nodes_explored < 1000) && - (time_since_log >= 1) || - (time_since_log > 60) || now > settings.time_limit) { - settings.log.printf(" %8d %8lu %+13.6e %+10.6e %4d %7.1e %s %9.2f\n", - nodes_explored, - heap.size(), - compute_user_objective(original_lp, upper_bound), - compute_user_objective(original_lp, lower_bound), - leaf_depth, - nodes_explored > 0 ? total_lp_iters / nodes_explored : 0, - user_mip_gap(compute_user_objective(original_lp, upper_bound), - compute_user_objective(original_lp, lower_bound)) - .c_str(), - now); - last_log = tic(); - } - if (now > settings.time_limit) { - settings.log.printf("Hit time limit. Stoppping\n"); - status = mip_status_t::TIME_LIMIT; - break; - } + mutex_branching_.lock(); + currently_branching_ = true; + mutex_branching_.unlock(); - // Set the correct bounds for the leaf problem - leaf_problem.lower = original_lp.lower; - leaf_problem.upper = original_lp.upper; - std::vector bounds_changed(original_lp.num_cols, false); - // Technically, we can get the already strengthened bounds from the node/parent instead of - // getting it from the original problem and re-strengthening. But this requires storing - // two vectors at each node and potentially cause memory issues - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - - std::vector& leaf_vstatus = node_ptr->vstatus; - lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); - - i_t node_iter = 0; - assert(leaf_vstatus.size() == leaf_problem.num_cols); - f_t lp_start_time = tic(); - std::vector leaf_edge_norms = edge_norms; // = node.steepest_edge_norms; - simplex_solver_settings_t lp_settings = settings; - lp_settings.set_log(false); - lp_settings.cut_off = upper_bound + settings.dual_tol; - lp_settings.inside_mip = 2; - - // in B&B we only have equality constraints, leave it empty for default - std::vector row_sense; - bool feasible = - bound_strengthening(row_sense, lp_settings, leaf_problem, Arow, var_types, bounds_changed); - - dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; - - // If the problem is infeasible after bounds strengthening, we don't need to solve the LP - if (feasible) { - lp_status = dual_phase2(2, - 0, - lp_start_time, - leaf_problem, - lp_settings, - leaf_vstatus, - leaf_solution, - node_iter, - leaf_edge_norms); - if (lp_status == dual::status_t::NUMERICAL) { - settings.log.printf("Numerical issue node %d. Resolving from scratch.\n", nodes_explored); - lp_status_t second_status = solve_linear_program_advanced( - leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); - lp_status = convert_lp_status_to_dual_status(second_status); - } - } - total_lp_solve_time += toc(lp_start_time); - total_lp_iters += node_iter; + std::future diving_thread; - nodes_explored++; - if (lp_status == dual::status_t::DUAL_UNBOUNDED) { - if (!feasible) { - settings.log.printf("Infeasible after bounds strengthening. Fathoming node %d.\n", - nodes_explored); - } - node_ptr->lower_bound = inf; - std::vector*> stack; - node_ptr->set_status(node_status_t::INFEASIBLE, stack); - graphviz_node(settings, node_ptr, "infeasible", 0.0); - remove_fathomed_nodes(stack); - // Node was infeasible. Do not branch - } else if (lp_status == dual::status_t::CUTOFF) { - node_ptr->lower_bound = upper_bound; - std::vector*> stack; - node_ptr->set_status(node_status_t::FATHOMED, stack); - f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); - graphviz_node(settings, node_ptr, "cut off", leaf_objective); - remove_fathomed_nodes(stack); - // Node was cut off. Do not branch - } else if (lp_status == dual::status_t::OPTIMAL) { - // LP was feasible - std::vector fractional; - const i_t leaf_fractional = - fractional_variables(settings, leaf_solution.x, var_types, fractional); - f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); - graphviz_node(settings, node_ptr, "lower bound", leaf_objective); - - pc.update_pseudo_costs(node_ptr, leaf_objective); - node_ptr->lower_bound = leaf_objective; - - constexpr f_t fathom_tol = 1e-5; - if (leaf_fractional == 0) { - bool send_solution = false; - mutex_upper.lock(); - if (leaf_objective < upper_bound_) { - incumbent.set_incumbent_solution(leaf_objective, leaf_solution.x); - upper_bound_ = upper_bound = leaf_objective; - gap = upper_bound - lower_bound; - settings.log.printf("B%8d %8lu %+13.6e %+10.6e %4d %7.1e %s %9.2f\n", - nodes_explored, - heap.size(), - compute_user_objective(original_lp, upper_bound), - compute_user_objective(original_lp, lower_bound), - leaf_depth, - nodes_explored > 0 ? total_lp_iters / nodes_explored : 0, - user_mip_gap(compute_user_objective(original_lp, upper_bound), - compute_user_objective(original_lp, lower_bound)) - .c_str(), - toc(start_time)); - send_solution = true; - } - mutex_upper.unlock(); - if (send_solution && settings.solution_callback != nullptr) { - std::vector original_x; - uncrush_primal_solution(original_problem, original_lp, incumbent.x, original_x); - settings.solution_callback(original_x, upper_bound); - } - graphviz_node(settings, node_ptr, "integer feasible", leaf_objective); - std::vector*> stack; - node_ptr->set_status(node_status_t::INTEGER_FEASIBLE, stack); - remove_fathomed_nodes(stack); - } else if (leaf_objective <= upper_bound + fathom_tol) { - // Choose fractional variable to branch on - const i_t branch_var = pc.variable_selection( - fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, log); - assert(leaf_vstatus.size() == leaf_problem.num_cols); - - // down child - std::unique_ptr> down_child = - std::make_unique>(original_lp, - node_ptr, - ++num_nodes, - branch_var, - 0, - leaf_solution.x[branch_var], - leaf_vstatus); - graphviz_edge(settings, - node_ptr, - down_child.get(), - branch_var, - 0, - std::floor(leaf_solution.x[branch_var])); - // up child - std::unique_ptr> up_child = - std::make_unique>(original_lp, - node_ptr, - ++num_nodes, - branch_var, - 1, - leaf_solution.x[branch_var], - leaf_vstatus); - graphviz_edge(settings, - node_ptr, - up_child.get(), - branch_var, - 0, - std::ceil(leaf_solution.x[branch_var])); - heap.push(down_child.get()); // the heap does not own the unique_ptr the tree does - heap.push(up_child.get()); // the heap does not own the unique_ptr the tree does - node_ptr->add_children(std::move(down_child), - std::move(up_child)); // child pointers moved into the tree - } else { - graphviz_node(settings, node_ptr, "fathomed", leaf_objective); - std::vector*> stack; - node_ptr->set_status(node_status_t::FATHOMED, stack); - remove_fathomed_nodes(stack); - } - } else { - graphviz_node(settings, node_ptr, "numerical", 0.0); - settings.log.printf("Encountered LP status %d. This indicates a numerical issue.\n", - lp_status); - status = mip_status_t::NUMERICAL; - break; - } + if (settings_.num_threads > 0) { + diving_thread = std::async(std::launch::async, [&]() { return dive(branch_var, solution); }); } - mutex_branching.lock(); - currently_branching = false; - mutex_branching.unlock(); - if (heap.size() == 0) { - mutex_lower.lock(); - lower_bound = lower_bound_ = root_node.lower_bound; - mutex_lower.unlock(); - gap = get_upper_bound() - lower_bound; - } + status = explore_tree(branch_var, solution); + + if (settings_.num_threads > 0) { mip_status_t diving_status = diving_thread.get(); } + + mutex_branching_.lock(); + currently_branching_ = false; + mutex_branching_.unlock(); - settings.log.printf( + settings_.log.printf( "Explored %d nodes in %.2fs.\nAbsolute Gap %e Objective %.16e Lower Bound %.16e\n", - nodes_explored, - toc(start_time), - gap, - compute_user_objective(original_lp, get_upper_bound()), - compute_user_objective(original_lp, lower_bound)); - - if (gap <= settings.absolute_mip_gap_tol || - relative_gap(get_upper_bound(), lower_bound) <= settings.relative_mip_gap_tol) { + stats_.nodes_explored.load(), + toc(stats_.start_time), + gap_, + compute_user_objective(original_lp_, get_upper_bound()), + compute_user_objective(original_lp_, lower_bound_)); + + if (gap_ <= settings_.absolute_mip_gap_tol || + relative_gap(get_upper_bound(), lower_bound_) <= settings_.relative_mip_gap_tol) { status = mip_status_t::OPTIMAL; - if (gap > 0 && gap <= settings.absolute_mip_gap_tol) { - settings.log.printf("Optimal solution found within absolute MIP gap tolerance (%.1e)\n", - settings.absolute_mip_gap_tol); - } else if (gap > 0 && - relative_gap(get_upper_bound(), lower_bound) <= settings.relative_mip_gap_tol) { - settings.log.printf("Optimal solution found within relative MIP gap tolerance (%.1e)\n", - settings.relative_mip_gap_tol); + if (gap_ > 0 && gap_ <= settings_.absolute_mip_gap_tol) { + settings_.log.printf("Optimal solution found within absolute MIP gap tolerance (%.1e)\n", + settings_.absolute_mip_gap_tol); + } else if (gap_ > 0 && + relative_gap(get_upper_bound(), lower_bound_) <= settings_.relative_mip_gap_tol) { + settings_.log.printf("Optimal solution found within relative MIP gap tolerance (%.1e)\n", + settings_.relative_mip_gap_tol); } else { - settings.log.printf("Optimal solution found.\n"); + settings_.log.printf("Optimal solution found.\n"); } - if (settings.heuristic_preemption_callback != nullptr) { - settings.heuristic_preemption_callback(); + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); } } - if (heap.size() == 0 && get_upper_bound() == inf) { - settings.log.printf("Integer infeasible.\n"); + if (stats_.nodes_unexplored == 0 && get_upper_bound() == inf) { + settings_.log.printf("Integer infeasible.\n"); status = mip_status_t::INFEASIBLE; - if (settings.heuristic_preemption_callback != nullptr) { - settings.heuristic_preemption_callback(); + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); } } - uncrush_primal_solution(original_problem, original_lp, incumbent.x, solution.x); - solution.objective = incumbent.objective; - solution.lower_bound = lower_bound; - solution.nodes_explored = nodes_explored; - solution.simplex_iterations = total_lp_iters; + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); + solution.objective = incumbent_.objective; + solution.lower_bound = get_lower_bound(); + solution.nodes_explored = stats_.nodes_explored; + solution.simplex_iterations = stats_.total_lp_iters; return status; } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 9c1bac9fb..efe28000b 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -18,15 +18,21 @@ #pragma once #include +#include #include +#include #include #include #include +#include +#include #include #include #include #include +#include "cuopt/linear_programming/mip/solver_settings.hpp" +#include "dual_simplex/mip_node.hpp" namespace cuopt::linear_programming::dual_simplex { @@ -50,65 +56,128 @@ class branch_and_bound_t { const simplex_solver_settings_t& solver_settings); // Set an initial guess based on the user_problem. This should be called before solve. - void set_initial_guess(const std::vector& user_guess) { guess = user_guess; } + void set_initial_guess(const std::vector& user_guess) { guess_ = user_guess; } // Set a solution based on the user problem during the course of the solve void set_new_solution(const std::vector& solution); - bool repair_solution(const std::vector& root_vstatus, - const std::vector& leaf_edge_norms, + // Repair a low-quality solution from the heuristics. + bool repair_solution(const std::vector& leaf_edge_norms, const std::vector& potential_solution, f_t& repaired_obj, std::vector& repaired_solution) const; f_t get_upper_bound(); + f_t get_lower_bound(); + // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); private: - const user_problem_t& original_problem; - const simplex_solver_settings_t settings; + const user_problem_t& original_problem_; + const simplex_solver_settings_t settings_; + + // Initial guess. + std::vector guess_; - f_t start_time; - std::vector guess; + lp_problem_t original_lp_; + std::vector new_slacks_; + std::vector var_types_; - lp_problem_t original_lp; - std::vector new_slacks; - std::vector var_types; // Mutex for lower bound - std::mutex mutex_lower; + std::mutex mutex_lower_; + // Global variable for lower bound f_t lower_bound_; // Mutex for upper bound - std::mutex mutex_upper; + std::mutex mutex_upper_; + // Global variable for upper bound f_t upper_bound_; + // Global variable for incumbent. The incumbent should be updated with the upper bound - mip_solution_t incumbent; + mip_solution_t incumbent_; // Mutex for gap - std::mutex mutex_gap; + std::mutex mutex_gap_; + // Global variable for gap - f_t gap; + f_t gap_; // Mutex for branching - std::mutex mutex_branching; - bool currently_branching; + std::mutex mutex_branching_; + bool currently_branching_; - // Mutex for stats - std::mutex mutex_stats; // Global variable for stats + std::mutex mutex_stats_; + + // Note that floating point atomics are only supported in C++20. struct stats_t { - int nodes_explored; - f_t total_lp_solve_time; - f_t start_time; - } stats; + f_t start_time = 0.0; + f_t total_lp_solve_time = 0.0; + std::atomic nodes_explored = 0; + std::atomic nodes_unexplored = 0; + f_t total_lp_iters = 0; + std::atomic num_nodes = 0; + } stats_; // Mutex for repair - std::mutex mutex_repair; - std::vector> repair_queue; + std::mutex mutex_repair_; + std::vector> repair_queue_; + + // Variables for the root node in the search tree. + std::vector root_vstatus_; + f_t root_objective_; + lp_solution_t root_relax_soln_; + std::vector edge_norms_; + + // Pseudocosts + pseudo_costs_t pc_; + std::mutex mutex_pc_; + + // Update the status of the nodes in the search tree. + void update_tree(mip_node_t* node_ptr, node_status_t status); + + // Update the incumbent solution with the new feasible solution. + // found during branch and bound. + void add_feasible_solution(f_t leaf_objective, + const std::vector& leaf_solution, + i_t leaf_depth, + char symbol); + + // Repairs low-quality solutions from the heuristics, if it is applicable. + void repair_heuristic_solutions(); + + // Explore the search tree using the best-first search strategy. + mip_status_t explore_tree(i_t branch_var, mip_solution_t& solution); + + // Explore the search tree using the depth-first search strategy. + mip_status_t dive(i_t branch_var, mip_solution_t& solution); + + // Branch the current node, creating two children. + void branch(mip_node_t* parent_node, + i_t branch_var, + f_t branch_var_val, + const std::vector& parent_vstatus); + + // Solve the LP relaxation of a leaf node. + mip_status_t solve_node_lp(mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + const std::vector& var_types, + f_t upper_bound); + + // Solve the LP relaxation of a leaf node using the dual simplex method. + dual::status_t node_dual_simplex(i_t leaf_id, + lp_problem_t& leaf_problem, + std::vector& leaf_vstatus, + lp_solution_t& leaf_solution, + std::vector& bounds_changed, + csc_matrix_t& Arow, + f_t upper_bound, + logger_t& log); }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index c4d83a833..2d315fe0e 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -33,6 +34,7 @@ enum class node_status_t : int { INTEGER_FEASIBLE = 2, // Node has an integer feasible solution INFEASIBLE = 3, // Node is infeasible FATHOMED = 4, // Node objective is greater than the upper bound + HAS_CHILDREN = 5, // Node has children to explore }; bool inactive_status(node_status_t status); @@ -40,7 +42,7 @@ bool inactive_status(node_status_t status); template class mip_node_t { public: - mip_node_t(f_t root_lower_bound, std::vector& basis) + mip_node_t(f_t root_lower_bound, const std::vector& basis) : status(node_status_t::ACTIVE), lower_bound(root_lower_bound), depth(0), @@ -53,13 +55,14 @@ class mip_node_t { children[0] = nullptr; children[1] = nullptr; } + mip_node_t(const lp_problem_t& problem, mip_node_t* parent_node, i_t node_num, i_t branch_variable, i_t branch_direction, f_t branch_var_value, - std::vector& basis) + const std::vector& basis) : status(node_status_t::ACTIVE), lower_bound(parent_node->lower_bound), depth(parent_node->depth + 1), @@ -103,6 +106,10 @@ class mip_node_t { } } + mip_node_t* get_down_child() const { return children[0].get(); } + + mip_node_t* get_up_child() const { return children[1].get(); } + void add_children(std::unique_ptr&& down_child, std::unique_ptr&& up_child) { diff --git a/cpp/src/dual_simplex/pseudo_costs.hpp b/cpp/src/dual_simplex/pseudo_costs.hpp index fef67b9b1..d26fe8d48 100644 --- a/cpp/src/dual_simplex/pseudo_costs.hpp +++ b/cpp/src/dual_simplex/pseudo_costs.hpp @@ -39,6 +39,14 @@ class pseudo_costs_t { void update_pseudo_costs(mip_node_t* node_ptr, f_t leaf_objective); + void resize(i_t num_variables) + { + pseudo_cost_sum_down.resize(num_variables); + pseudo_cost_sum_up.resize(num_variables); + pseudo_cost_num_down.resize(num_variables); + pseudo_cost_num_up.resize(num_variables); + } + void initialized(i_t& num_initialized_down, i_t& num_initialized_up, f_t& pseudo_cost_down_avg, diff --git a/cpp/src/mip/diversity/recombiners/sub_mip.cuh b/cpp/src/mip/diversity/recombiners/sub_mip.cuh index eb9435351..1818092fc 100644 --- a/cpp/src/mip/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip/diversity/recombiners/sub_mip.cuh @@ -17,6 +17,7 @@ #pragma once +#include "cuopt/linear_programming/mip/solver_settings.hpp" #include "recombiner.cuh" #include @@ -115,6 +116,7 @@ class sub_mip_recombiner_t : public recombiner_t { f_t objective) { this->solution_callback(solution, objective); }; + // disable B&B logs, so that it is not interfering with the main B&B thread branch_and_bound_settings.log.log = false; dual_simplex::branch_and_bound_t branch_and_bound(branch_and_bound_problem, diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp index 923ad8c90..964ea314c 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve.cpp +++ b/cpp/tests/dual_simplex/unit_tests/solve.cpp @@ -159,6 +159,7 @@ TEST(dual_simplex, burglar) for (int j = 0; j < num_items; ++j) { user_problem.var_types[j] = cuopt::linear_programming::dual_simplex::variable_type_t::INTEGER; } + cuopt::linear_programming::dual_simplex::simplex_solver_settings_t settings; std::vector solution(num_items); EXPECT_EQ((cuopt::linear_programming::dual_simplex::solve(user_problem, settings, solution)), 0); diff --git a/cpp/tests/mip/miplib_test.cu b/cpp/tests/mip/miplib_test.cu index 603fe6014..d0866455f 100644 --- a/cpp/tests/mip/miplib_test.cu +++ b/cpp/tests/mip/miplib_test.cu @@ -16,6 +16,9 @@ */ #include "../linear_programming/utilities/pdlp_test_utilities.cuh" +#include "cuopt/linear_programming/mip/solver_settings.hpp" +#include "dual_simplex/branch_and_bound.hpp" +#include "dual_simplex/simplex_solver_settings.hpp" #include "mip_utils.cuh" #include @@ -40,7 +43,7 @@ struct result_map_t { double cost; }; -void test_miplib_file(result_map_t test_instance) +void test_miplib_file(result_map_t test_instance, mip_solver_settings_t settings) { const raft::handle_t handle_{}; @@ -48,7 +51,6 @@ void test_miplib_file(result_map_t test_instance) cuopt::mps_parser::mps_data_model_t problem = cuopt::mps_parser::parse_mps(path, false); handle_.sync_stream(); - mip_solver_settings_t settings; // set the time limit depending on we are in assert mode or not #ifdef ASSERT_MODE constexpr double test_time_limit = 60.; @@ -58,7 +60,9 @@ void test_miplib_file(result_map_t test_instance) settings.time_limit = test_time_limit; mip_solution_t solution = solve_mip(&handle_, problem, settings); - EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::FeasibleFound); + bool is_feasible = solution.get_termination_status() == mip_termination_status_t::FeasibleFound || + solution.get_termination_status() == mip_termination_status_t::Optimal; + EXPECT_TRUE(is_feasible); double obj_val = solution.get_objective_value(); // for now keep a 100% error rate EXPECT_NEAR(test_instance.cost, obj_val, test_instance.cost); @@ -68,10 +72,11 @@ void test_miplib_file(result_map_t test_instance) TEST(mip_solve, run_small_tests) { + mip_solver_settings_t settings; std::vector test_instances = { {"mip/50v-10.mps", 11311031.}, {"mip/neos5.mps", 15.}, {"mip/swath1.mps", 1300.}}; for (const auto& test_instance : test_instances) { - test_miplib_file(test_instance); + test_miplib_file(test_instance, settings); } }