diff --git a/benchmarks/linear_programming/cuopt/mip_test_instances.hpp b/benchmarks/linear_programming/cuopt/mip_test_instances.hpp index 5ca15cbb41..ac7c673590 100644 --- a/benchmarks/linear_programming/cuopt/mip_test_instances.hpp +++ b/benchmarks/linear_programming/cuopt/mip_test_instances.hpp @@ -17,66 +17,175 @@ #pragma once #include #include -std::vector instances = { - // "30n20b8.mps", - // "bab2.mps", - // "bab6.mps", - // "bnatt400.mps", - // "bnatt500.mps", // infeasible - // "cmflsp50-24-8-8.mps", - // "cryptanalysiskb128n5obj14.mps", //infeasible - // "cryptanalysiskb128n5obj16.mps", - // "csched007.mps", - // "csched008.mps", - "dano3_3.mps", - // "dano3_5.mps", - // "dws008-01.mps", - // "enlight_hard.mps", - // "fhnw-binpack4-48.mps", - // "fhnw-binpack4-4.mps", // infeasible - "gfd-schedulen180f7d50m30k18.mps", - // "highschool1-aigio.mps", - "icir97_tension.mps", - "irish-electricity.mps", - // "istanbul-no-cutoff.mps", - // "milo-v12-6-r2-40-1.mps", - // "momentum1.mps", - // "neos-1354092.mps", - // "neos-2075418-temuka.mps", // infeasible - // "neos-3024952-loue.mps", - // "neos-3216931-puriri.mps", - // "neos-3381206-awhea.mps", - // "neos-3402454-bohle.mps", // infeasible - // "neos-3656078-kumeu.mps", - // "neos-3988577-wolgan.mps", // infeasible - // "neos-4413714-turia.mps", - "neos-4647030-tutaki.mps", - "neos-4338804-snowy.mps", - "neos-4722843-widden.mps", - "neos-4763324-toguru.mps", - "neos-5104907-jarama.mps", - "neos-5114902-kasavu.mps", - // "neos859080.mps", // infeasible - "ns1760995.mps", - "ns1952667.mps", - // "nursesched-medium-hint03.mps", - // "nursesched-sprint02.mps", - // "peg-solitaire-a3.mps", - // "physiciansched3-3.mps", - // "physiciansched6-2.mps", - "piperout-08.mps", - // "piperout-27.mps", - // "rail01.mps", - // "rail02.mps", - // "rd-rplusc-21.mps", - // "rocI-4-11.mps", - "rocII-5-11.mps", - // "roll3000.mps", - "s100.mps", - // "supportcase19.mps", - // "supportcase22.mps", - "swath1.mps" - // "swath3.mps", - // "timtab1.mps", - // "triptim1.mps" -}; +std::vector instances = {"30n20b8.mps", + "50v-10.mps", + "CMS750_4.mps", + "academictimetablesmall.mps", + "air05.mps", + "app1-1.mps", + "app1-2.mps", + "assign1-5-8.mps", + "atlanta-ip.mps", + "bab2.mps", + "bab6.mps", + "beasleyC3.mps", + "binkar10_1.mps", + "blp-ar98.mps", + "blp-ic98.mps", + "bppc4-08.mps", + "brazil3.mps", + "cmflsp50-24-8-8.mps", + "co-100.mps", + "cod105.mps", + "comp07-2idx.mps", + "comp21-2idx.mps", + "csched007.mps", + "csched008.mps", + "cvs16r128-89.mps", + "dano3_3.mps", + "decomp2.mps", + "drayage-100-23.mps", + "drayage-25-23.mps", + "eil33-2.mps", + "eilA101-2.mps", + "exp-1-500-5-5.mps", + "fast0507.mps", + "fastxgemm-n2r6s0t2.mps", + "fiball.mps", + "gen-ip002.mps", + "germanrr.mps", + "glass4.mps", + "graph20-20-1rand.mps", + "graphdraw-domain.mps", + "h80x6320d.mps", + "highschool1-aigio.mps", + "hypothyroid-k1.mps", + "icir97_tension.mps", + "irish-electricity.mps", + "istanbul-no-cutoff.mps", + "k1mushroom.mps", + "lectsched-5-obj.mps", + "leo1.mps", + "leo2.mps", + "lotsize.mps", + "mad.mps", + "map10.mps", + "map16715-04.mps", + "markshare2.mps", + "markshare_4_0.mps", + "mas74.mps", + "mc11.mps", + "mcsched.mps", + "mik-250-20-75-4.mps", + "momentum1.mps", + "mushroom-best.mps", + "mzzv11.mps", + "mzzv42z.mps", + "n2seq36q.mps", + "n3div36.mps", + "neos-1171448.mps", + "neos-1171737.mps", + "neos-1354092.mps", + "neos-1445765.mps", + "neos-1456979.mps", + "neos-1582420.mps", + "neos-2657525-crna.mps", + "neos-2746589-doon.mps", + "neos-3024952-loue.mps", + "neos-3046615-murg.mps", + "neos-3216931-puriri.mps", + "neos-3402294-bobin.mps", + "neos-3656078-kumeu.mps", + "neos-3754480-nidda.mps", + "neos-4300652-rahue.mps", + "neos-4338804-snowy.mps", + "neos-4387871-tavua.mps", + "neos-4413714-turia.mps", + "neos-4532248-waihi.mps", + "neos-4722843-widden.mps", + "neos-4738912-atrato.mps", + "neos-4763324-toguru.mps", + "neos-4954672-berkel.mps", + "neos-5049753-cuanza.mps", + "neos-5093327-huahum.mps", + "neos-5107597-kakapo.mps", + "neos-5114902-kasavu.mps", + "neos-5188808-nattai.mps", + "neos-5195221-niemur.mps", + "neos-662469.mps", + "neos-787933.mps", + "neos-848589.mps", + "neos-860300.mps", + "neos-911970.mps", + "neos-933966.mps", + "neos-950242.mps", + "neos17.mps", + "neos5.mps", + "net12.mps", + "netdiversion.mps", + "nexp-150-20-8-5.mps", + "ns1644855.mps", + "ns1760995.mps", + "ns1830653.mps", + "nursesched-medium-hint03.mps", + "nursesched-sprint02.mps", + "opm2-z10-s4.mps", + "pg.mps", + "physiciansched3-3.mps", + "piperout-08.mps", + "piperout-27.mps", + "pk1.mps", + "qap10.mps", + "radiationm18-12-05.mps", + "radiationm40-10-02.mps", + "rail01.mps", + "rail02.mps", + "rail507.mps", + "ran14x18-disj-8.mps", + "rmatr100-p10.mps", + "rmatr200-p5.mps", + "rocI-4-11.mps", + "rocII-5-11.mps", + "rococoB10-011000.mps", + "rococoC10-001000.mps", + "roi2alpha3n4.mps", + "roi5alpha10n8.mps", + "roll3000.mps", + "s100.mps", + "s250r10.mps", + "satellites2-40.mps", + "satellites2-60-fs.mps", + "savsched1.mps", + "sct2.mps", + "seymour.mps", + "seymour1.mps", + "sing326.mps", + "sing44.mps", + "sorrell3.mps", + "sp97ar.mps", + "sp98ar.mps", + "splice1k1.mps", + "square41.mps", + "square47.mps", + "supportcase10.mps", + "supportcase12.mps", + "supportcase18.mps", + "supportcase26.mps", + "supportcase33.mps", + "supportcase40.mps", + "supportcase42.mps", + "supportcase6.mps", + "supportcase7.mps", + "swath1.mps", + "swath3.mps", + "tbfp-network.mps", + "thor50dday.mps", + "timtab1.mps", + "tr12-30.mps", + "traininstance2.mps", + "traininstance6.mps", + "trento1.mps", + "uccase12.mps", + "uct-subprob.mps", + "unitcal_7.mps", + "var-smallemery-m6j6.mps"}; diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp index a2243f5405..55eb0f656c 100644 --- a/cpp/src/dual_simplex/basis_solves.cpp +++ b/cpp/src/dual_simplex/basis_solves.cpp @@ -349,7 +349,6 @@ i_t factorize_basis(const csc_matrix_t& A, } } } - assert(Snz <= Snz_max && (Sdim > 0 && Snz > 0)); S.col_start[Sdim] = Snz; // Finalize S csc_matrix_t SL(Sdim, Sdim, Snz); diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index ecd7488f19..ce92532044 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -30,55 +30,11 @@ #include #include -#include -#include #include #include namespace cuopt::linear_programming::dual_simplex { -namespace global_variables { - -#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE - -// Mutex for lower bound -std::mutex mutex_lower; -// Global variable for lower bound -double lower_bound; - -// Mutex for upper bound -std::mutex mutex_upper; -// Global variable for upper bound -double upper_bound; -// Global variable for incumbent. The incumbent should be updated with the upper bound -mip_solution_t incumbent(1); - -// Mutex for gap -std::mutex mutex_gap; -// Global variable for gap -double gap; - -// Mutex for branching -std::mutex mutex_branching; -bool currently_branching; - -// Mutex for stats -std::mutex mutex_stats; -// Global variable for stats -struct stats_t { - int nodes_explored; - double total_lp_solve_time; - double start_time; -} stats; - -// Mutex for repair -std::mutex mutex_repair; -std::vector> repair_queue; - -#endif - -} // namespace global_variables - namespace { template @@ -230,12 +186,12 @@ dual::status_t convert_lp_status_to_dual_status(lp_status_t status) } // namespace -template -f_t get_upper_bound() +template +f_t branch_and_bound_t::get_upper_bound() { - global_variables::mutex_upper.lock(); - const f_t upper_bound = global_variables::upper_bound; - global_variables::mutex_upper.unlock(); + mutex_upper.lock(); + const f_t upper_bound = upper_bound_; + mutex_upper.unlock(); return upper_bound; } @@ -282,15 +238,15 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu f_t obj = compute_objective(original_lp, crushed_solution); bool is_feasible = false; bool attempt_repair = false; - global_variables::mutex_upper.lock(); - if (obj < global_variables::upper_bound) { + 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); if (is_feasible) { - global_variables::upper_bound = obj; + upper_bound_ = obj; } else { attempt_repair = true; constexpr bool verbose = false; @@ -304,15 +260,15 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu } } } - global_variables::mutex_upper.unlock(); + mutex_upper.unlock(); if (is_feasible) { - global_variables::mutex_lower.lock(); - f_t lower_bound = global_variables::lower_bound; - global_variables::mutex_lower.unlock(); - global_variables::mutex_branching.lock(); - bool currently_branching = global_variables::currently_branching; - global_variables::mutex_branching.unlock(); + mutex_lower.lock(); + f_t lower_bound = lower_bound_; + mutex_lower.unlock(); + mutex_branching.lock(); + bool currently_branching = currently_branching; + mutex_branching.unlock(); if (currently_branching) { settings.log.printf( "H %+13.6e %+10.6e %s %9.2f\n", @@ -330,9 +286,9 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu } if (attempt_repair) { - global_variables::mutex_repair.lock(); - global_variables::repair_queue.push_back(crushed_solution); - global_variables::mutex_repair.unlock(); + mutex_repair.lock(); + repair_queue.push_back(crushed_solution); + mutex_repair.unlock(); } } @@ -400,32 +356,28 @@ 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) + : original_problem(user_problem), settings(solver_settings), original_lp(1, 1, 1), incumbent(1) { start_time = tic(); convert_user_problem(original_problem, settings, original_lp, new_slacks); full_variable_types(original_problem, original_lp, var_types); - global_variables::mutex_upper.lock(); - global_variables::upper_bound = inf; - global_variables::mutex_upper.unlock(); + mutex_upper.lock(); + upper_bound_ = inf; + mutex_upper.unlock(); - global_variables::mutex_lower.lock(); - global_variables::lower_bound = -inf; - global_variables::mutex_lower.unlock(); + mutex_lower.lock(); + lower_bound_ = -inf; + mutex_lower.unlock(); - global_variables::mutex_branching.lock(); - global_variables::currently_branching = false; - global_variables::mutex_branching.unlock(); + mutex_branching.lock(); + currently_branching = false; + mutex_branching.unlock(); } template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { - // TODO remove this after submip PR merge - global_variables::mutex_upper.lock(); - global_variables::repair_queue.clear(); - global_variables::mutex_upper.unlock(); mip_status_t status = mip_status_t::UNSET; mip_solution_t incumbent(original_lp.num_cols); if (guess.size() != 0) { @@ -438,10 +390,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut 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); - global_variables::mutex_upper.lock(); + mutex_upper.lock(); incumbent.set_incumbent_solution(computed_obj, crushed_guess); - global_variables::upper_bound = computed_obj; - global_variables::mutex_upper.unlock(); + upper_bound_ = computed_obj; + mutex_upper.unlock(); } } @@ -479,21 +431,21 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut 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.solution_callback != nullptr) { + if (settings.set_simplex_solution_callback != nullptr) { std::vector original_x; uncrush_primal_solution(original_problem, original_lp, root_relax_soln.x, original_x); settings.set_simplex_solution_callback(original_x, compute_user_objective(original_lp, root_objective)); } - global_variables::mutex_lower.lock(); - f_t lower_bound = global_variables::lower_bound = root_objective; - global_variables::mutex_lower.unlock(); + mutex_lower.lock(); + f_t lower_bound = lower_bound_ = root_objective; + mutex_lower.unlock(); if (num_fractional == 0) { - global_variables::mutex_upper.lock(); + mutex_upper.lock(); incumbent.set_incumbent_solution(root_objective, root_relax_soln.x); - global_variables::upper_bound = root_objective; - global_variables::mutex_upper.unlock(); + 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; @@ -577,28 +529,28 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut 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 - f_t gap = get_upper_bound() - lower_bound; + f_t gap = get_upper_bound() - lower_bound; i_t nodes_explored = 0; settings.log.printf( "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap | " " Time \n"); - global_variables::mutex_branching.lock(); - global_variables::currently_branching = true; - global_variables::mutex_branching.unlock(); + 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 && + 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; - global_variables::mutex_repair.lock(); - if (global_variables::repair_queue.size() > 0) { - to_repair = global_variables::repair_queue; - global_variables::repair_queue.clear(); + mutex_repair.lock(); + if (repair_queue.size() > 0) { + to_repair = repair_queue; + repair_queue.clear(); } - global_variables::mutex_repair.unlock(); + mutex_repair.unlock(); if (to_repair.size() > 0) { settings.log.debug("Attempting to repair %ld injected solutions\n", to_repair.size()); @@ -608,9 +560,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut bool is_feasible = repair_solution( root_vstatus, edge_norms, potential_solution, repaired_obj, repaired_solution); if (is_feasible) { - global_variables::mutex_upper.lock(); - if (repaired_obj < global_variables::upper_bound) { - global_variables::upper_bound = repaired_obj; + mutex_upper.lock(); + if (repaired_obj < upper_bound_) { + upper_bound_ = repaired_obj; incumbent.set_incumbent_solution(repaired_obj, repaired_solution); settings.log.printf( @@ -627,7 +579,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings.solution_callback(original_x, repaired_obj); } } - global_variables::mutex_upper.unlock(); + mutex_upper.unlock(); } } } @@ -635,7 +587,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // 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(); + 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 @@ -645,9 +597,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut remove_fathomed_nodes(stack); continue; } - global_variables::mutex_lower.lock(); - global_variables::lower_bound = lower_bound = node_ptr->lower_bound; - global_variables::mutex_lower.unlock(); + 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); @@ -740,11 +692,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut constexpr f_t fathom_tol = 1e-5; if (leaf_fractional == 0) { bool send_solution = false; - global_variables::mutex_upper.lock(); - if (leaf_objective < global_variables::upper_bound) { + mutex_upper.lock(); + if (leaf_objective < upper_bound_) { incumbent.set_incumbent_solution(leaf_objective, leaf_solution.x); - global_variables::upper_bound = upper_bound = leaf_objective; - gap = upper_bound - lower_bound; + 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(), @@ -758,7 +710,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut toc(start_time)); send_solution = true; } - global_variables::mutex_upper.unlock(); + 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); @@ -822,15 +774,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut break; } } - global_variables::mutex_branching.lock(); - global_variables::currently_branching = false; - global_variables::mutex_branching.unlock(); + mutex_branching.lock(); + currently_branching = false; + mutex_branching.unlock(); if (heap.size() == 0) { - global_variables::mutex_lower.lock(); - lower_bound = global_variables::lower_bound = root_node.lower_bound; - global_variables::mutex_lower.unlock(); - gap = get_upper_bound() - lower_bound; + mutex_lower.lock(); + lower_bound = lower_bound_ = root_node.lower_bound; + mutex_lower.unlock(); + gap = get_upper_bound() - lower_bound; } settings.log.printf( @@ -838,17 +790,17 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut nodes_explored, toc(start_time), gap, - compute_user_objective(original_lp, get_upper_bound()), + 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) { + 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) { + 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 { @@ -859,7 +811,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } } - if (heap.size() == 0 && get_upper_bound() == inf) { + if (heap.size() == 0 && get_upper_bound() == inf) { settings.log.printf("Integer infeasible.\n"); status = mip_status_t::INFEASIBLE; if (settings.heuristic_preemption_callback != nullptr) { diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index da594bf93b..9c1bac9fbc 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -23,6 +23,8 @@ #include #include +#include +#include #include #include @@ -59,6 +61,8 @@ class branch_and_bound_t { f_t& repaired_obj, std::vector& repaired_solution) const; + f_t get_upper_bound(); + // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); @@ -72,6 +76,39 @@ class branch_and_bound_t { lp_problem_t original_lp; std::vector new_slacks; std::vector var_types; + // Mutex for lower bound + std::mutex mutex_lower; + // Global variable for lower bound + f_t lower_bound_; + + // Mutex for upper bound + 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; + + // Mutex for gap + std::mutex mutex_gap; + // Global variable for gap + f_t gap; + + // Mutex for branching + std::mutex mutex_branching; + bool currently_branching; + + // Mutex for stats + std::mutex mutex_stats; + // Global variable for stats + struct stats_t { + int nodes_explored; + f_t total_lp_solve_time; + f_t start_time; + } stats; + + // Mutex for repair + std::mutex mutex_repair; + std::vector> repair_queue; }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/pseudo_costs.cpp b/cpp/src/dual_simplex/pseudo_costs.cpp index d351d7f704..95c96f859f 100644 --- a/cpp/src/dual_simplex/pseudo_costs.cpp +++ b/cpp/src/dual_simplex/pseudo_costs.cpp @@ -21,27 +21,10 @@ #include #include -#include #include namespace cuopt::linear_programming::dual_simplex { -namespace global_variables { - -#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE - -// Global variable storing the results of strong branching down -std::vector strong_branch_down; -// Global variable storing the results of strong branching up -std::vector strong_branch_up; - -std::mutex strong_branches_completed; -int num_strong_branches_completed = 0; - -#endif - -} // namespace global_variables - namespace { template @@ -56,7 +39,8 @@ void strong_branch_helper(i_t thread_id, f_t root_obj, const std::vector& root_soln, const std::vector& root_vstatus, - const std::vector& edge_norms) + const std::vector& edge_norms, + pseudo_costs_t& pc) { lp_problem_t child_problem = original_lp; @@ -112,7 +96,7 @@ void strong_branch_helper(i_t thread_id, } if (branch == 0) { - global_variables::strong_branch_down[k] = obj - root_obj; + pc.strong_branch_down[k] = obj - root_obj; if (verbose) { settings.log.printf("Thread id %2d remaining %d variable %d branch %d obj %e time %.2f\n", thread_id, @@ -123,7 +107,7 @@ void strong_branch_helper(i_t thread_id, toc(start_time)); } } else { - global_variables::strong_branch_up[k] = obj - root_obj; + pc.strong_branch_up[k] = obj - root_obj; if (verbose) { settings.log.printf( "Thread id %2d remaining %d variable %d branch %d obj %e change down %e change up %e " @@ -133,8 +117,8 @@ void strong_branch_helper(i_t thread_id, j, branch, obj, - global_variables::strong_branch_down[k], - global_variables::strong_branch_up[k], + pc.strong_branch_down[k], + pc.strong_branch_up[k], toc(start_time)); } } @@ -142,15 +126,15 @@ void strong_branch_helper(i_t thread_id, } if (toc(start_time) > settings.time_limit) { break; } - global_variables::strong_branches_completed.lock(); - global_variables::num_strong_branches_completed++; - global_variables::strong_branches_completed.unlock(); + pc.strong_branches_completed.lock(); + pc.num_strong_branches_completed++; + pc.strong_branches_completed.unlock(); if (thread_id == 0 && toc(last_log) > 10) { last_log = tic(); - global_variables::strong_branches_completed.lock(); - const i_t completed = global_variables::num_strong_branches_completed; - global_variables::strong_branches_completed.unlock(); + pc.strong_branches_completed.lock(); + const i_t completed = pc.num_strong_branches_completed; + pc.strong_branches_completed.unlock(); settings.log.printf("%d of %ld strong branches completed in %.1fs\n", completed, fractional.size(), @@ -190,7 +174,8 @@ void thread_strong_branch(i_t thread_id, f_t root_obj, const std::vector& root_soln, const std::vector& root_vstatus, - const std::vector& edge_norms) + const std::vector& edge_norms, + pseudo_costs_t& pc) { i_t start, end; get_partitioning(static_cast(fractional.size()), num_threads, thread_id, start, end); @@ -211,7 +196,8 @@ void thread_strong_branch(i_t thread_id, root_obj, root_soln, root_vstatus, - edge_norms); + edge_norms, + pc); } } // namespace @@ -228,10 +214,10 @@ void strong_branching(const lp_problem_t& original_lp, const std::vector& edge_norms, pseudo_costs_t& pc) { - global_variables::strong_branch_down.resize(fractional.size()); - global_variables::strong_branch_up.resize(fractional.size()); + pc.strong_branch_down.resize(fractional.size()); + pc.strong_branch_up.resize(fractional.size()); - global_variables::num_strong_branches_completed = 0; + pc.num_strong_branches_completed = 0; if (fractional.size() > settings.num_threads) { int num_threads = settings.num_threads; @@ -251,7 +237,8 @@ void strong_branching(const lp_problem_t& original_lp, root_obj, root_soln, root_vstatus, - edge_norms); + edge_norms, + std::ref(pc)); } for (int thread_id = 0; thread_id < num_threads; thread_id++) { @@ -270,13 +257,11 @@ void strong_branching(const lp_problem_t& original_lp, root_obj, root_soln, root_vstatus, - edge_norms); + edge_norms, + pc); } - pc.update_pseudo_costs_from_strong_branching(fractional, - root_soln, - global_variables::strong_branch_down, - global_variables::strong_branch_up); + pc.update_pseudo_costs_from_strong_branching(fractional, root_soln); } template @@ -393,10 +378,7 @@ i_t pseudo_costs_t::variable_selection(const std::vector& fractio template void pseudo_costs_t::update_pseudo_costs_from_strong_branching( - const std::vector& fractional, - const std::vector& root_soln, - const std::vector& strong_branch_down, - const std::vector& strong_branch_up) + const std::vector& fractional, const std::vector& root_soln) { for (i_t k = 0; k < fractional.size(); k++) { const i_t j = fractional[k]; diff --git a/cpp/src/dual_simplex/pseudo_costs.hpp b/cpp/src/dual_simplex/pseudo_costs.hpp index 74be1ce2b7..fef67b9b14 100644 --- a/cpp/src/dual_simplex/pseudo_costs.hpp +++ b/cpp/src/dual_simplex/pseudo_costs.hpp @@ -22,6 +22,8 @@ #include #include +#include + namespace cuopt::linear_programming::dual_simplex { template @@ -49,13 +51,16 @@ class pseudo_costs_t { logger_t& log) const; void update_pseudo_costs_from_strong_branching(const std::vector& fractional, - const std::vector& root_soln, - const std::vector& strong_branch_down, - const std::vector& strong_branch_up); + const std::vector& root_soln); std::vector pseudo_cost_sum_up; std::vector pseudo_cost_sum_down; std::vector pseudo_cost_num_up; std::vector pseudo_cost_num_down; + std::vector strong_branch_down; + std::vector strong_branch_up; + + std::mutex strong_branches_completed; + i_t num_strong_branches_completed = 0; }; template diff --git a/cpp/src/mip/diversity/diversity_config.hpp b/cpp/src/mip/diversity/diversity_config.hpp index abdbdf0bb8..c38555ab90 100644 --- a/cpp/src/mip/diversity/diversity_config.hpp +++ b/cpp/src/mip/diversity/diversity_config.hpp @@ -20,26 +20,25 @@ namespace cuopt::linear_programming::detail { struct diversity_config_t { - static constexpr double time_ratio_on_init_lp = 0.1; - static constexpr double max_time_on_lp = 30; - static constexpr double time_ratio_of_probing_cache = 0.10; - static constexpr double max_time_on_probing = 60; - static constexpr size_t max_iterations_without_improvement = 15; - static constexpr int max_var_diff = 256; - static constexpr size_t max_solutions = 32; - static constexpr double initial_infeasibility_weight = 1000.; - static constexpr double default_time_limit = 10.; - static constexpr int initial_island_size = 3; - static constexpr int maximum_island_size = 8; - static constexpr bool use_avg_diversity = false; - static constexpr double generation_time_limit_ratio = 0.6; - static constexpr double max_island_gen_time = 600; - static constexpr size_t n_sol_for_skip_init_gen = 3; - static constexpr double max_fast_sol_time = 10; - static constexpr double lp_run_time_if_feasible = 15.; - static constexpr double lp_run_time_if_infeasible = 1; - static constexpr double close_to_parents_ratio = 0.1; - static constexpr bool halve_population = false; + double time_ratio_on_init_lp = 0.1; + double max_time_on_lp = 30; + double time_ratio_of_probing_cache = 0.10; + double max_time_on_probing = 60; + size_t max_iterations_without_improvement = 15; + int max_var_diff = 256; + size_t max_solutions = 32; + double initial_infeasibility_weight = 1000.; + double default_time_limit = 10.; + int initial_island_size = 3; + int maximum_island_size = 8; + bool use_avg_diversity = false; + double generation_time_limit_ratio = 0.6; + double max_island_gen_time = 600; + size_t n_sol_for_skip_init_gen = 3; + double max_fast_sol_time = 10; + double lp_run_time_if_feasible = 15.; + double lp_run_time_if_infeasible = 1; + bool halve_population = true; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu index e0d96d447c..b406d56a39 100644 --- a/cpp/src/mip/diversity/diversity_manager.cu +++ b/cpp/src/mip/diversity/diversity_manager.cu @@ -19,7 +19,6 @@ #include #include #include -#include "diversity_config.hpp" #include "diversity_manager.cuh" #include @@ -32,34 +31,29 @@ constexpr bool fp_only_run = false; namespace cuopt::linear_programming::detail { -constexpr int max_var_diff = diversity_config_t::max_var_diff; -constexpr size_t max_solutions = diversity_config_t::max_solutions; -constexpr double initial_infeasibility_weight = diversity_config_t::initial_infeasibility_weight; -constexpr double default_time_limit = diversity_config_t::default_time_limit; -constexpr int initial_island_size = diversity_config_t::initial_island_size; -constexpr int maximum_island_size = diversity_config_t::maximum_island_size; -constexpr bool use_avg_diversity = diversity_config_t::use_avg_diversity; - size_t fp_recombiner_config_t::max_n_of_vars_from_other = fp_recombiner_config_t::initial_n_of_vars_from_other; size_t ls_recombiner_config_t::max_n_of_vars_from_other = ls_recombiner_config_t::initial_n_of_vars_from_other; size_t bp_recombiner_config_t::max_n_of_vars_from_other = bp_recombiner_config_t::initial_n_of_vars_from_other; +size_t sub_mip_recombiner_config_t::max_n_of_vars_from_other = + sub_mip_recombiner_config_t::initial_n_of_vars_from_other; template diversity_manager_t::diversity_manager_t(mip_solver_context_t& context_) - : problem_ptr(context.problem_ptr), - context(context_), + : context(context_), + problem_ptr(context.problem_ptr), + diversity_config(), population("population", context, - max_var_diff, - max_solutions, - initial_infeasibility_weight * context.problem_ptr->n_constraints), + diversity_config.max_var_diff, + diversity_config.max_solutions, + diversity_config.initial_infeasibility_weight * context.problem_ptr->n_constraints), lp_optimal_solution(context.problem_ptr->n_variables, context.problem_ptr->handle_ptr->get_stream()), ls(context, lp_optimal_solution), - timer(default_time_limit), + timer(diversity_config.default_time_limit), bound_prop_recombiner(context, context.problem_ptr->n_variables, ls.constraint_prop, @@ -70,10 +64,14 @@ diversity_manager_t::diversity_manager_t(mip_solver_context_tn_variables, ls.line_segment_search, context.problem_ptr->handle_ptr), + sub_mip_recombiner( + context, population, context.problem_ptr->n_variables, context.problem_ptr->handle_ptr), rng(cuopt::seed_generator::get_seed()), stats(context.stats), - mab_recombiner( - recombiner_enum_t::SIZE, cuopt::seed_generator::get_seed(), recombiner_alpha, "recombiner"), + mab_recombiner(static_cast(recombiner_enum_t::SIZE), + cuopt::seed_generator::get_seed(), + recombiner_alpha, + "recombiner"), mab_ls(mab_ls_config_t::n_of_arms, cuopt::seed_generator::get_seed(), ls_alpha, "ls"), assignment_hash_map(*context.problem_ptr) { @@ -84,7 +82,7 @@ diversity_manager_t::diversity_manager_t(mip_solver_context_t::diversity_manager_t(mip_solver_context_t::generate_initial_solutions() { add_user_given_solutions(initial_sol_vector); bool skip_initial_island_generation = - initial_sol_vector.size() > diversity_config_t::n_sol_for_skip_init_gen || from_dir; + initial_sol_vector.size() > diversity_config.n_sol_for_skip_init_gen || from_dir; // allocate maximum of 40% of the time to the initial island generation // aim to generate at least 5 feasible solutions thus spending 8% of the time to generate a // solution if we can generate faster generate up to 10 sols const f_t generation_time_limit = - diversity_config_t::generation_time_limit_ratio * timer.get_time_limit(); - const f_t max_island_gen_time = diversity_config_t::max_island_gen_time; + diversity_config.generation_time_limit_ratio * timer.get_time_limit(); + const f_t max_island_gen_time = diversity_config.max_island_gen_time; f_t total_island_gen_time = std::min(generation_time_limit, max_island_gen_time); timer_t gen_timer(total_island_gen_time); f_t sol_time_limit = gen_timer.remaining_time(); - for (i_t i = 0; i < maximum_island_size && !skip_initial_island_generation; ++i) { + for (i_t i = 0; i < diversity_config.maximum_island_size && !skip_initial_island_generation; + ++i) { if (check_b_b_preemption()) { return; } if (i + population.get_external_solution_size() >= 5) { break; } CUOPT_LOG_DEBUG("Generating sol %d", i); bool is_first_sol = (i == 0); - if (i == 1) { sol_time_limit = gen_timer.remaining_time() / (initial_island_size - 1); } + if (i == 1) { + sol_time_limit = gen_timer.remaining_time() / (diversity_config.initial_island_size - 1); + } // in first iteration, definitely generate feasible if (is_first_sol) { sol_time_limit = gen_timer.remaining_time(); @@ -280,13 +278,13 @@ void diversity_manager_t::generate_initial_solutions() average_fj_weights(i); // run ls on the solutions // if at least initial_island_size solutions are generated and time limit is reached - if (i >= initial_island_size || gen_timer.check_time_limit()) { break; } + if (i >= diversity_config.initial_island_size || gen_timer.check_time_limit()) { break; } } CUOPT_LOG_DEBUG("Initial unsearched solutions are generated!"); i_t actual_island_size = initial_sol_vector.size(); population.normalize_weights(); // find diversity of the population - population.find_diversity(initial_sol_vector, use_avg_diversity); + population.find_diversity(initial_sol_vector, diversity_config.use_avg_diversity); population.add_solutions_from_vec(std::move(initial_sol_vector)); population.update_qualities(); CUOPT_LOG_DEBUG("Initial population generated, size %d var_threshold %d", @@ -332,7 +330,7 @@ void diversity_manager_t::generate_quick_feasible_solution() solution_t solution(*problem_ptr); // min 1 second, max 10 seconds const f_t generate_fast_solution_time = - std::min(diversity_config_t::max_fast_sol_time, std::max(1., timer.remaining_time() / 20.)); + std::min(diversity_config.max_fast_sol_time, std::max(1., timer.remaining_time() / 20.)); timer_t sol_timer(generate_fast_solution_time); // do very short LP run to get somewhere close to the optimal point ls.generate_fast_solution(solution, sol_timer); @@ -395,10 +393,10 @@ void diversity_manager_t::run_fp_alone(solution_t& solution) template solution_t diversity_manager_t::run_solver() { - population.timer = timer; - const f_t time_limit = timer.remaining_time(); - const f_t lp_time_limit = std::min(diversity_config_t::max_time_on_lp, - time_limit * diversity_config_t::time_ratio_on_init_lp); + population.timer = timer; + const f_t time_limit = timer.remaining_time(); + const f_t lp_time_limit = + std::min(diversity_config.max_time_on_lp, time_limit * diversity_config.time_ratio_on_init_lp); // to automatically compute the solving time on scope exit auto timer_raii_guard = cuopt::scope_guard([&]() { stats.total_solve_time = timer.elapsed_time(); }); @@ -421,8 +419,8 @@ solution_t diversity_manager_t::run_solver() if (check_b_b_preemption()) { return population.best_feasible(); } // before probing cache or LP, run FJ to generate initial primal feasible solution if (!from_dir && !fp_only_run && !fj_only_run) { generate_quick_feasible_solution(); } - constexpr f_t time_ratio_of_probing_cache = diversity_config_t::time_ratio_of_probing_cache; - constexpr f_t max_time_on_probing = diversity_config_t::max_time_on_probing; + const f_t time_ratio_of_probing_cache = diversity_config.time_ratio_of_probing_cache; + const f_t max_time_on_probing = diversity_config.max_time_on_probing; f_t time_for_probing_cache = std::min(max_time_on_probing, time_limit * time_ratio_of_probing_cache); timer_t probing_timer{time_for_probing_cache}; @@ -515,8 +513,8 @@ template void diversity_manager_t::diversity_step() { // TODO when the solver is faster, increase this number - constexpr i_t max_iterations_without_improvement = - diversity_config_t::max_iterations_without_improvement; + const i_t max_iterations_without_improvement = + diversity_config.max_iterations_without_improvement; bool improved = true; while (improved) { int k = max_iterations_without_improvement; @@ -625,7 +623,7 @@ void diversity_manager_t::main_loop() diversity_step(); if (timer.check_time_limit()) { break; } - if (diversity_config_t::halve_population) { + if (diversity_config.halve_population) { population.adjust_threshold(timer); i_t prev_threshold = population.var_threshold; population.halve_the_population(); @@ -635,7 +633,7 @@ void diversity_manager_t::main_loop() current_population.insert(current_population.end(), std::make_move_iterator(new_solutions.begin()), std::make_move_iterator(new_solutions.end())); - population.find_diversity(current_population, use_avg_diversity); + population.find_diversity(current_population, diversity_config.use_avg_diversity); // if the threshold is lower than the threshold we progress with time // set it to the higher threshold // population.var_threshold = max(population.var_threshold, prev_threshold); @@ -694,7 +692,7 @@ diversity_manager_t::recombine_and_local_search(solution_t& auto [offspring, success] = recombine(sol1, sol2); if (!success) { // add the attempt - mab_recombiner.add_mab_reward(recombine_stats.get_last_attempt(), + mab_recombiner.add_mab_reward(static_cast(recombine_stats.get_last_attempt()), std::numeric_limits::lowest(), std::numeric_limits::lowest(), std::numeric_limits::max(), @@ -714,7 +712,7 @@ diversity_manager_t::recombine_and_local_search(solution_t& success = this->run_local_search(offspring, population.weights, timer, ls_config); if (!success) { // add the attempt - mab_recombiner.add_mab_reward(recombine_stats.get_last_attempt(), + mab_recombiner.add_mab_reward(static_cast(recombine_stats.get_last_attempt()), std::numeric_limits::lowest(), std::numeric_limits::lowest(), std::numeric_limits::max(), @@ -732,8 +730,8 @@ diversity_manager_t::recombine_and_local_search(solution_t& solution_t lp_offspring(offspring); cuopt_assert(population.test_invariant(), ""); cuopt_assert(lp_offspring.test_number_all_integer(), "All must be integers before LP"); - f_t lp_run_time = offspring.get_feasible() ? diversity_config_t::lp_run_time_if_feasible - : diversity_config_t::lp_run_time_if_infeasible; + f_t lp_run_time = offspring.get_feasible() ? diversity_config.lp_run_time_if_feasible + : diversity_config.lp_run_time_if_infeasible; lp_run_time = std::min(lp_run_time, timer.remaining_time()); relaxed_lp_settings_t lp_settings; lp_settings.time_limit = lp_run_time; @@ -758,16 +756,16 @@ diversity_manager_t::recombine_and_local_search(solution_t& f_t best_quality_of_parents = std::min(sol1.get_quality(population.weights), sol2.get_quality(population.weights)); mab_recombiner.add_mab_reward( - recombine_stats.get_last_attempt(), + static_cast(recombine_stats.get_last_attempt()), best_quality_of_parents, population.best().get_quality(population.weights), offspring_qual, recombiner_work_normalized_reward_t(recombine_stats.get_last_recombiner_time())); - mab_ls.add_mab_reward(mab_ls_config_t::last_ls_mab_option, - best_quality_of_parents, - population.best_feasible().get_quality(population.weights), - offspring_qual, - ls_work_normalized_reward_t(mab_ls_config_t::last_lm_config)); + // mab_ls.add_mab_reward(mab_ls_config_t::last_ls_mab_option, + // best_quality_of_parents, + // population.best_feasible().get_quality(population.weights), + // offspring_qual, + // ls_work_normalized_reward_t(mab_ls_config_t::last_lm_config)); if (context.settings.benchmark_info_ptr != nullptr) { check_better_than_both(offspring, sol1, sol2); check_better_than_both(lp_offspring, sol1, sol2); @@ -779,34 +777,53 @@ template std::pair, bool> diversity_manager_t::recombine( solution_t& a, solution_t& b) { - i_t recombiner; + recombiner_enum_t recombiner; if (run_only_ls_recombiner) { recombiner = recombiner_enum_t::LINE_SEGMENT; } else if (run_only_bp_recombiner) { recombiner = recombiner_enum_t::BOUND_PROP; } else if (run_only_fp_recombiner) { recombiner = recombiner_enum_t::FP; + } else if (run_only_sub_mip_recombiner) { + recombiner = recombiner_enum_t::SUB_MIP; } else { - recombiner = mab_recombiner.select_mab_option(); + recombiner = static_cast(mab_recombiner.select_mab_option()); } recombine_stats.add_attempt((recombiner_enum_t)recombiner); recombine_stats.start_recombiner_time(); - if (recombiner == recombiner_enum_t::BOUND_PROP) { - auto [sol, success] = bound_prop_recombiner.recombine(a, b, population.weights); - recombine_stats.stop_recombiner_time(); - if (success) { recombine_stats.add_success(); } - return std::make_pair(sol, success); - } else if (recombiner == recombiner_enum_t::FP) { - auto [sol, success] = fp_recombiner.recombine(a, b, population.weights); - recombine_stats.stop_recombiner_time(); - if (success) { recombine_stats.add_success(); } - return std::make_pair(sol, success); - } else { - auto [sol, success] = line_segment_recombiner.recombine(a, b, population.weights); - recombine_stats.stop_recombiner_time(); - if (success) { recombine_stats.add_success(); } - return std::make_pair(sol, success); + // Refactored code using a switch statement + switch (recombiner) { + case recombiner_enum_t::BOUND_PROP: { + auto [sol, success] = bound_prop_recombiner.recombine(a, b, population.weights); + recombine_stats.stop_recombiner_time(); + if (success) { recombine_stats.add_success(); } + return std::make_pair(sol, success); + } + case recombiner_enum_t::FP: { + auto [sol, success] = fp_recombiner.recombine(a, b, population.weights); + recombine_stats.stop_recombiner_time(); + if (success) { recombine_stats.add_success(); } + return std::make_pair(sol, success); + } + case recombiner_enum_t::LINE_SEGMENT: { + auto [sol, success] = line_segment_recombiner.recombine(a, b, population.weights); + recombine_stats.stop_recombiner_time(); + if (success) { recombine_stats.add_success(); } + return std::make_pair(sol, success); + } + case recombiner_enum_t::SUB_MIP: { + auto [sol, success] = sub_mip_recombiner.recombine(a, b, population.weights); + recombine_stats.stop_recombiner_time(); + if (success) { recombine_stats.add_success(); } + return std::make_pair(sol, success); + } + case recombiner_enum_t::SIZE: { + CUOPT_LOG_ERROR("Invalid or unhandled recombiner type: %d", recombiner); + return std::make_pair(solution_t(a), false); + } } + CUOPT_LOG_ERROR("Invalid or unhandled recombiner type: %d", recombiner); + return std::make_pair(solution_t(a), false); } template diff --git a/cpp/src/mip/diversity/diversity_manager.cuh b/cpp/src/mip/diversity/diversity_manager.cuh index caf4d68036..3ad538d33b 100644 --- a/cpp/src/mip/diversity/diversity_manager.cuh +++ b/cpp/src/mip/diversity/diversity_manager.cuh @@ -18,6 +18,7 @@ #pragma once #include "assignment_hash_map.cuh" +#include "diversity_config.hpp" #include "multi_armed_bandit.cuh" #include "population.cuh" @@ -25,6 +26,7 @@ #include "recombiners/fp_recombiner.cuh" #include "recombiners/line_segment_recombiner.cuh" #include "recombiners/recombiner_stats.hpp" +#include "recombiners/sub_mip.cuh" #include #include @@ -79,6 +81,7 @@ class diversity_manager_t { mip_solver_context_t& context; problem_t* problem_ptr; + diversity_config_t diversity_config; population_t population; rmm::device_uvector lp_optimal_solution; bool simplex_solution_exists{false}; @@ -87,6 +90,7 @@ class diversity_manager_t { bound_prop_recombiner_t bound_prop_recombiner; fp_recombiner_t fp_recombiner; line_segment_recombiner_t line_segment_recombiner; + sub_mip_recombiner_t sub_mip_recombiner; all_recombine_stats recombine_stats; std::mt19937 rng; i_t current_step{0}; @@ -103,6 +107,7 @@ class diversity_manager_t { bool run_only_ls_recombiner{false}; bool run_only_bp_recombiner{false}; bool run_only_fp_recombiner{false}; + bool run_only_sub_mip_recombiner{false}; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/diversity/multi_armed_bandit.cuh b/cpp/src/mip/diversity/multi_armed_bandit.cuh index d520f3202f..d4783069f5 100644 --- a/cpp/src/mip/diversity/multi_armed_bandit.cuh +++ b/cpp/src/mip/diversity/multi_armed_bandit.cuh @@ -30,7 +30,6 @@ constexpr double ls_alpha = 0.03; template struct mab_ls_config_t { - // RANDOM is the last ls method static constexpr i_t n_of_arms = static_cast(ls_method_t::LS_METHODS_SIZE); static void get_local_search_and_lm_from_config(i_t config_id, ls_config_t& ls_config) diff --git a/cpp/src/mip/diversity/recombiners/recombiner_configs.hpp b/cpp/src/mip/diversity/recombiners/recombiner_configs.hpp index 18c9590150..59ac5c3c24 100644 --- a/cpp/src/mip/diversity/recombiners/recombiner_configs.hpp +++ b/cpp/src/mip/diversity/recombiners/recombiner_configs.hpp @@ -33,14 +33,16 @@ struct bp_recombiner_config_t { static void increase_max_n_of_vars_from_other() { max_n_of_vars_from_other = std::min( - size_t(max_n_of_vars_from_other * n_var_ratio_increase_factor), max_different_var_limit); + static_cast(std::ceil(max_n_of_vars_from_other * n_var_ratio_increase_factor)), + max_different_var_limit); CUOPT_LOG_DEBUG("Increased max_n_of_vars_from_other in BP recombiner to %lu", max_n_of_vars_from_other); } static void decrease_max_n_of_vars_from_other() { max_n_of_vars_from_other = std::max( - size_t(max_n_of_vars_from_other * n_var_ratio_decrease_factor), min_different_var_limit); + static_cast(std::floor(max_n_of_vars_from_other * n_var_ratio_decrease_factor)), + min_different_var_limit); CUOPT_LOG_DEBUG("Decreased max_n_of_vars_from_other in BP recombiner to %lu", max_n_of_vars_from_other); } @@ -61,14 +63,16 @@ struct ls_recombiner_config_t { static void increase_max_n_of_vars_from_other() { max_n_of_vars_from_other = std::min( - size_t(max_n_of_vars_from_other * n_var_ratio_increase_factor), max_different_var_limit); + static_cast(std::ceil(max_n_of_vars_from_other * n_var_ratio_increase_factor)), + max_different_var_limit); CUOPT_LOG_DEBUG("Increased max_n_of_vars_from_other in LS recombiner to %lu", max_n_of_vars_from_other); } static void decrease_max_n_of_vars_from_other() { max_n_of_vars_from_other = std::max( - size_t(max_n_of_vars_from_other * n_var_ratio_decrease_factor), min_different_var_limit); + static_cast(std::floor(max_n_of_vars_from_other * n_var_ratio_decrease_factor)), + min_different_var_limit); CUOPT_LOG_DEBUG("Decreased max_n_of_vars_from_other in LS recombiner to %lu", max_n_of_vars_from_other); } @@ -88,17 +92,46 @@ struct fp_recombiner_config_t { static void increase_max_n_of_vars_from_other() { max_n_of_vars_from_other = std::min( - size_t(max_n_of_vars_from_other * n_var_ratio_increase_factor), max_different_var_limit); + static_cast(std::ceil(max_n_of_vars_from_other * n_var_ratio_increase_factor)), + max_different_var_limit); CUOPT_LOG_DEBUG("Increased max_n_of_vars_from_other in FP recombiner to %lu", max_n_of_vars_from_other); } static void decrease_max_n_of_vars_from_other() { max_n_of_vars_from_other = std::max( - size_t(max_n_of_vars_from_other * n_var_ratio_decrease_factor), min_different_var_limit); + static_cast(std::floor(max_n_of_vars_from_other * n_var_ratio_decrease_factor)), + min_different_var_limit); CUOPT_LOG_DEBUG("Decreased max_n_of_vars_from_other in FP recombiner to %lu", max_n_of_vars_from_other); } }; +struct sub_mip_recombiner_config_t { + static constexpr double sub_mip_time_limit = 3.; + static constexpr double infeasibility_detection_time_limit = 0.5; + static constexpr size_t initial_n_of_vars_from_other = 40; + static constexpr size_t max_different_var_limit = 500; + static constexpr size_t min_different_var_limit = 10; + static size_t max_n_of_vars_from_other; + static constexpr double n_var_ratio_increase_factor = 1.1; + static constexpr double n_var_ratio_decrease_factor = 0.99; + static void increase_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::min( + static_cast(std::ceil(max_n_of_vars_from_other * n_var_ratio_increase_factor)), + max_different_var_limit); + CUOPT_LOG_DEBUG("Increased max_n_of_vars_from_other in SUB_MIP recombiner to %lu", + max_n_of_vars_from_other); + } + static void decrease_max_n_of_vars_from_other() + { + max_n_of_vars_from_other = std::max( + static_cast(std::floor(max_n_of_vars_from_other * n_var_ratio_decrease_factor)), + min_different_var_limit); + CUOPT_LOG_DEBUG("Decreased max_n_of_vars_from_other in SUB_MIP recombiner to %lu", + max_n_of_vars_from_other); + } +}; + } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp b/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp index 1486999880..7e3e25164a 100644 --- a/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp +++ b/cpp/src/mip/diversity/recombiners/recombiner_stats.hpp @@ -21,7 +21,7 @@ namespace cuopt::linear_programming::detail { -enum recombiner_enum_t : int { BOUND_PROP = 0, FP, LINE_SEGMENT, SIZE }; +enum class recombiner_enum_t : int { BOUND_PROP = 0, FP, LINE_SEGMENT, SUB_MIP, SIZE }; struct recombine_stats { int attempts; @@ -65,7 +65,7 @@ struct recombine_stats { struct all_recombine_stats { static constexpr size_t recombiner_count = static_cast(recombiner_enum_t::SIZE); - static constexpr std::array recombiner_labels = {"BOUND_PROP", "FP", "LINE_SEGMENT"}; + static constexpr std::array recombiner_labels = {"BOUND_PROP", "FP", "LINE_SEGMENT", "SUB_MIP"}; std::array stats; diff --git a/cpp/src/mip/diversity/recombiners/sub_mip.cuh b/cpp/src/mip/diversity/recombiners/sub_mip.cuh new file mode 100644 index 0000000000..381f5b1258 --- /dev/null +++ b/cpp/src/mip/diversity/recombiners/sub_mip.cuh @@ -0,0 +1,185 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#include "recombiner.cuh" + +#include +#include +#include + +namespace cuopt::linear_programming::detail { + +template +class sub_mip_recombiner_t : public recombiner_t { + public: + sub_mip_recombiner_t(mip_solver_context_t& context, + population_t& population, + i_t n_vars, + const raft::handle_t* handle_ptr) + : recombiner_t(context, n_vars, handle_ptr), + vars_to_fix(n_vars, handle_ptr->get_stream()), + context(context), + population(population) + { + } + + void solution_callback(std::vector& solution, f_t objective) + { + CUOPT_LOG_DEBUG("SUBMIP added solution with objective %.16e", objective); + solution_vector.push_back(solution); + } + + std::pair, bool> recombine(solution_t& a, + solution_t& b, + const weight_t& weights) + { + raft::common::nvtx::range fun_scope("Sub-MIP recombiner"); + solution_vector.clear(); + auto& guiding_solution = a.get_feasible() ? a : b; + auto& other_solution = a.get_feasible() ? b : a; + // copy the solution from A + solution_t offspring(guiding_solution); + // find same values and populate it to offspring + i_t n_different_vars = + this->assign_same_integer_values(guiding_solution, other_solution, offspring); + CUOPT_LOG_DEBUG("SUB_MIP rec: Number of different variables %d MAX_VARS %d", + n_different_vars, + sub_mip_recombiner_config_t::max_n_of_vars_from_other); + i_t n_vars_from_other = n_different_vars; + if (n_vars_from_other > (i_t)sub_mip_recombiner_config_t::max_n_of_vars_from_other) { + n_vars_from_other = sub_mip_recombiner_config_t::max_n_of_vars_from_other; + thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()}; + thrust::shuffle(a.handle_ptr->get_thrust_policy(), + this->remaining_indices.data(), + this->remaining_indices.data() + n_different_vars, + g); + } + i_t n_vars_from_guiding = a.problem_ptr->n_integer_vars - n_vars_from_other; + if (n_vars_from_other == 0 || n_vars_from_guiding == 0) { + CUOPT_LOG_DEBUG("Returning false because all vars are common or different"); + return std::make_pair(offspring, false); + } + CUOPT_LOG_DEBUG( + "n_vars_from_guiding %d n_vars_from_other %d", n_vars_from_guiding, n_vars_from_other); + this->compute_vars_to_fix(offspring, vars_to_fix, n_vars_from_other, n_vars_from_guiding); + auto [fixed_problem, fixed_assignment, variable_map] = offspring.fix_variables(vars_to_fix); + fixed_problem.presolve_data.reset_additional_vars(fixed_problem, offspring.handle_ptr); + fixed_problem.presolve_data.initialize_var_mapping(fixed_problem, offspring.handle_ptr); + trivial_presolve(fixed_problem); + fixed_problem.check_problem_representation(true); + // brute force rounding threshold is 8 + const bool run_sub_mip = fixed_problem.n_integer_vars > 8; + dual_simplex::mip_status_t branch_and_bound_status = dual_simplex::mip_status_t::UNSET; + dual_simplex::mip_solution_t branch_and_bound_solution(1); + if (run_sub_mip) { + // run sub-mip + namespace dual_simplex = cuopt::linear_programming::dual_simplex; + dual_simplex::user_problem_t branch_and_bound_problem; + dual_simplex::simplex_solver_settings_t branch_and_bound_settings; + fixed_problem.get_host_user_problem(branch_and_bound_problem); + branch_and_bound_solution.resize(branch_and_bound_problem.num_cols); + // Fill in the settings for branch and bound + branch_and_bound_settings.time_limit = sub_mip_recombiner_config_t::sub_mip_time_limit; + branch_and_bound_settings.print_presolve_stats = false; + branch_and_bound_settings.absolute_mip_gap_tol = context.settings.tolerances.absolute_mip_gap; + branch_and_bound_settings.relative_mip_gap_tol = context.settings.tolerances.relative_mip_gap; + branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; + branch_and_bound_settings.solution_callback = [this](std::vector& solution, + 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, + branch_and_bound_settings); + branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution); + if (solution_vector.size() > 0) { + cuopt_assert(fixed_assignment.size() == branch_and_bound_solution.x.size(), + "Assignment size mismatch"); + CUOPT_LOG_DEBUG("Sub-MIP solution found. Objective %.16e. Status %d", + branch_and_bound_solution.objective, + int(branch_and_bound_status)); + // first post process the trivial presolve on a device vector + rmm::device_uvector post_processed_solution(branch_and_bound_solution.x.size(), + offspring.handle_ptr->get_stream()); + raft::copy(post_processed_solution.data(), + branch_and_bound_solution.x.data(), + branch_and_bound_solution.x.size(), + offspring.handle_ptr->get_stream()); + fixed_problem.post_process_assignment(post_processed_solution, false); + cuopt_assert(post_processed_solution.size() == fixed_assignment.size(), + "Assignment size mismatch"); + offspring.handle_ptr->sync_stream(); + std::swap(fixed_assignment, post_processed_solution); + } + offspring.handle_ptr->sync_stream(); + } + if (solution_vector.size() > 0) { + // unfix the assignment on given result no matter if it is feasible + offspring.unfix_variables(fixed_assignment, variable_map); + } else { + offspring.round_nearest(); + } + cuopt_func_call(offspring.test_variable_bounds()); + cuopt_assert(offspring.test_number_all_integer(), "All must be integers after offspring"); + offspring.compute_feasibility(); + // bool same_as_parents = this->check_if_offspring_is_same_as_parents(offspring, a, b); + // adjust the max_n_of_vars_from_other + if (n_different_vars > (i_t)sub_mip_recombiner_config_t::max_n_of_vars_from_other) { + if (branch_and_bound_status == dual_simplex::mip_status_t::OPTIMAL) { + sub_mip_recombiner_config_t::increase_max_n_of_vars_from_other(); + } else { + sub_mip_recombiner_config_t::decrease_max_n_of_vars_from_other(); + } + } + // try adding all intermediate solutions to the population + for (const auto& solution : solution_vector) { + CUOPT_LOG_DEBUG("Adding intermediate submip solution to population"); + solution_t sol(offspring); + rmm::device_uvector fixed_assignment(solution.size(), + offspring.handle_ptr->get_stream()); + raft::copy(fixed_assignment.data(), + solution.data(), + solution.size(), + offspring.handle_ptr->get_stream()); + fixed_problem.post_process_assignment(fixed_assignment, false); + sol.unfix_variables(fixed_assignment, variable_map); + sol.compute_feasibility(); + cuopt_func_call(sol.test_variable_bounds()); + population.add_solution(std::move(sol)); + } + bool better_cost_than_parents = + offspring.get_quality(weights) < + std::min(other_solution.get_quality(weights), guiding_solution.get_quality(weights)); + bool better_feasibility_than_parents = offspring.get_feasible() && + !other_solution.get_feasible() && + !guiding_solution.get_feasible(); + if (better_cost_than_parents || better_feasibility_than_parents) { + CUOPT_LOG_DEBUG("Offspring is feasible or better than both parents"); + return std::make_pair(offspring, true); + } + return std::make_pair(offspring, !std::isnan(branch_and_bound_solution.objective)); + } + rmm::device_uvector vars_to_fix; + mip_solver_context_t& context; + std::vector> solution_vector; + population_t& population; +}; + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/local_search/local_search.cu b/cpp/src/mip/local_search/local_search.cu index 5e39491264..ae8a416fa9 100644 --- a/cpp/src/mip/local_search/local_search.cu +++ b/cpp/src/mip/local_search/local_search.cu @@ -125,6 +125,7 @@ bool local_search_t::run_local_search(solution_t& solution, } else { fj.copy_weights(weights, solution.handle_ptr); is_feas = run_fj_annealing(solution, timer, ls_config); + if (lp_optimal_exists) { is_feas = run_fj_line_segment(solution, timer, ls_config); } } return is_feas; } diff --git a/cpp/src/mip/local_search/rounding/constraint_prop.cu b/cpp/src/mip/local_search/rounding/constraint_prop.cu index e557d3ba81..6e6a5deb37 100644 --- a/cpp/src/mip/local_search/rounding/constraint_prop.cu +++ b/cpp/src/mip/local_search/rounding/constraint_prop.cu @@ -599,11 +599,10 @@ thrust::pair constraint_prop_t::generate_double_probing_pair f_t random_value = dist(rng); f_t average_value = (from_first + from_second) / 2; if (random_value > 0.5) { - average_value = ceil(average_value); + average_value = ceil(average_value - orig_sol.problem_ptr->tolerances.integrality_tolerance); } else { - average_value = floor(average_value); + average_value = floor(average_value + orig_sol.problem_ptr->tolerances.integrality_tolerance); } - // if the first set of values are more represented use the second value if (probing_config.value().get().n_of_fixed_from_first > probing_config.value().get().n_of_fixed_from_second) { diff --git a/cpp/src/mip/problem/presolve_data.cuh b/cpp/src/mip/problem/presolve_data.cuh index 218d5f077e..20ae5fe931 100644 --- a/cpp/src/mip/problem/presolve_data.cuh +++ b/cpp/src/mip/problem/presolve_data.cuh @@ -19,11 +19,15 @@ #include +#include #include namespace cuopt { namespace linear_programming::detail { +template +class problem_t; + template class presolve_data_t { public: @@ -37,6 +41,7 @@ class presolve_data_t { fixed_var_assignment(0, stream) { } + presolve_data_t(const presolve_data_t& other, rmm::cuda_stream_view stream) : variable_offsets(other.variable_offsets), additional_var_used(other.additional_var_used), @@ -47,6 +52,26 @@ class presolve_data_t { fixed_var_assignment(other.fixed_var_assignment, stream) { } + + void initialize_var_mapping(const problem_t& problem, const raft::handle_t* handle_ptr) + { + variable_mapping.resize(problem.n_variables, handle_ptr->get_stream()); + thrust::sequence( + handle_ptr->get_thrust_policy(), variable_mapping.begin(), variable_mapping.end()); + fixed_var_assignment.resize(problem.n_variables, handle_ptr->get_stream()); + thrust::uninitialized_fill(handle_ptr->get_thrust_policy(), + fixed_var_assignment.begin(), + fixed_var_assignment.end(), + 0.); + } + + void reset_additional_vars(const problem_t& problem, const raft::handle_t* handle_ptr) + { + variable_offsets.assign(problem.n_variables, 0); + additional_var_used.assign(problem.n_variables, false); + additional_var_id_per_var.assign(problem.n_variables, -1); + } + presolve_data_t(presolve_data_t&&) = default; presolve_data_t& operator=(presolve_data_t&&) = default; presolve_data_t& operator=(const presolve_data_t&) = delete; diff --git a/cpp/src/mip/problem/problem.cu b/cpp/src/mip/problem/problem.cu index b3ca958d2d..1a5f76b038 100644 --- a/cpp/src/mip/problem/problem.cu +++ b/cpp/src/mip/problem/problem.cu @@ -86,8 +86,6 @@ void problem_t::op_problem_cstr_body(const optimization_problem_tget_stream()); is_binary_variable.resize(n_variables, handle_ptr->get_stream()); compute_n_integer_vars(); @@ -670,7 +668,8 @@ bool problem_t::pre_process_assignment(rmm::device_uvector& assig // it removes the additional variable for free variables // and expands the assignment to the original variable dimension template -void problem_t::post_process_assignment(rmm::device_uvector& current_assignment) +void problem_t::post_process_assignment(rmm::device_uvector& current_assignment, + bool resize_to_original_problem) { cuopt_assert(current_assignment.size() == presolve_data.variable_mapping.size(), "size mismatch"); auto assgn = make_span(current_assignment); @@ -700,7 +699,9 @@ void problem_t::post_process_assignment(rmm::device_uvector& curr raft::copy( current_assignment.data(), h_assignment.data(), h_assignment.size(), handle_ptr->get_stream()); // this separate resizing is needed because of the callback - current_assignment.resize(original_problem_ptr->get_n_variables(), handle_ptr->get_stream()); + if (resize_to_original_problem) { + current_assignment.resize(original_problem_ptr->get_n_variables(), handle_ptr->get_stream()); + } } template @@ -1470,15 +1471,7 @@ void problem_t::preprocess_problem() compute_csr(variable_constraint_map, *this); compute_transpose_of_problem(); check_problem_representation(true, false); - presolve_data.variable_mapping.resize(n_variables, handle_ptr->get_stream()); - thrust::sequence(handle_ptr->get_thrust_policy(), - presolve_data.variable_mapping.begin(), - presolve_data.variable_mapping.end()); - presolve_data.fixed_var_assignment.resize(n_variables, handle_ptr->get_stream()); - thrust::uninitialized_fill(handle_ptr->get_thrust_policy(), - presolve_data.fixed_var_assignment.begin(), - presolve_data.fixed_var_assignment.end(), - 0.); + presolve_data.initialize_var_mapping(*this, handle_ptr); integer_indices.resize(n_variables, handle_ptr->get_stream()); is_binary_variable.resize(n_variables, handle_ptr->get_stream()); original_ids.resize(n_variables); diff --git a/cpp/src/mip/problem/problem.cuh b/cpp/src/mip/problem/problem.cuh index d8b57ac6f3..49103fc955 100644 --- a/cpp/src/mip/problem/problem.cuh +++ b/cpp/src/mip/problem/problem.cuh @@ -95,7 +95,8 @@ class problem_t { void resize_constraints(size_t matrix_size, size_t constraint_size, size_t var_size); void preprocess_problem(); bool pre_process_assignment(rmm::device_uvector& assignment); - void post_process_assignment(rmm::device_uvector& current_assignment); + void post_process_assignment(rmm::device_uvector& current_assignment, + bool resize_to_original_problem = true); void post_process_solution(solution_t& solution); void compute_transpose_of_problem(); f_t get_user_obj_from_solver_obj(f_t solver_obj);