From 930447c76c230581dad3960d1142dd3677905ca6 Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Tue, 13 May 2025 14:34:49 +0200 Subject: [PATCH 1/2] set_cover: backport from main --- ortools/set_cover/samples/cft.cc | 53 ++ ortools/set_cover/set_cover_cft.cc | 732 +++++++++++++++++------- ortools/set_cover/set_cover_cft.h | 132 +++-- ortools/set_cover/set_cover_submodel.cc | 105 ++-- ortools/set_cover/set_cover_submodel.h | 35 +- ortools/set_cover/set_cover_views.h | 10 + ortools/set_cover/views.h | 7 + 7 files changed, 778 insertions(+), 296 deletions(-) create mode 100644 ortools/set_cover/samples/cft.cc diff --git a/ortools/set_cover/samples/cft.cc b/ortools/set_cover/samples/cft.cc new file mode 100644 index 0000000000..90a72ea6b0 --- /dev/null +++ b/ortools/set_cover/samples/cft.cc @@ -0,0 +1,53 @@ +#include +#include +#include + +#include "ortools/base/init_google.h" +#include "ortools/set_cover/set_cover_cft.h" +#include "ortools/set_cover/set_cover_reader.h" + +using namespace operations_research; +ABSL_FLAG(std::string, instance, "", "SCP instance int RAIL format."); + +#define DO_PRICING +int main(int argc, char **argv) { + InitGoogle(argv[0], &argc, &argv, true); + + scp::Model original_model = ReadOrlibRail(absl::GetFlag(FLAGS_instance)); + +#ifdef DO_PRICING + scp::FullToCoreModel model(&original_model); +#else + scp::SubModel model(&original_model); +#endif + + scp::PrimalDualState result = scp::RunCftHeuristic(model); + auto [solution, dual] = result; + if (solution.subsets().empty()) { + std::cerr << "Error: failed to find any solution\n"; + } else { + std::cout << "Solution: " << solution.cost() << '\n'; + } + +#ifdef DO_PRICING + if (dual.multipliers().empty()) { + std::cerr << "Error: failed to find any dual\n"; + } else { + std::cout << "Core Lower bound: " << dual.lower_bound() << '\n'; + } + if (model.best_dual_state().multipliers().empty()) { + std::cerr << "Error: no real dual state has been computed\n"; + } else { + std::cout << "Full Lower bound: " << model.best_dual_state().lower_bound() + << '\n'; + } +#else + if (dual.multipliers().empty()) { + std::cerr << "Error: failed to find any dual\n"; + } else { + std::cout << "Lower bound: " << dual.lower_bound() << '\n'; + } +#endif + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 71628ed1b3..88af8d42d5 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -14,39 +14,73 @@ #include "ortools/set_cover/set_cover_cft.h" #include +#include +#include #include #include #include +#include #include +#include #include #include -#include "ortools/base/status_macros.h" #include "ortools/base/stl_util.h" #include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_submodel.h" #include "ortools/set_cover/set_cover_views.h" +#include "ortools/set_cover/views.h" namespace operations_research::scp { +// Minimum distance between lower and upper bounds to consider them different. +// If cost are all integral, can be set neear to 1.0 +#define CFT_BOUND_EPSILON .999 +#define CFT_MAX_MULTIPLIER 1e9 +#define CFT_MEASURE_TIME + //////////////////////////////////////////////////////////////////////// ////////////////////////// COMMON DEFINITIONS ////////////////////////// //////////////////////////////////////////////////////////////////////// namespace { -class CoverCounters { + +#ifdef CFT_MEASURE_TIME +#define CFT_MEASURE_SCOPE_DURATION(Timer) \ + Timer.Start(); \ + Defer pause_timer = [&] { Timer.Stop(); }; + +thread_local WallTimer subgradient_time; +thread_local WallTimer greedy_time; +thread_local WallTimer three_phase_time; +thread_local WallTimer refinement_time; +#else +#define CFT_MEASURE_SCOPE_DURATION(Timer) +#endif + +// StopWatch does not add up duration of multiple invocations, Defer is a lower +// level construct useful in this case. +template +struct Defer : T { + Defer(T&& t) : T(std::forward(t)) {} + ~Defer() { T::operator()(); } +}; + +template +class CoverCountersImpl { public: - CoverCounters(BaseInt nelems = 0) : cov_counters(nelems, 0) {} + CoverCountersImpl(BaseInt nelems = 0) : cov_counters(nelems, 0) {} void Reset(BaseInt nelems) { cov_counters.assign(nelems, 0); } BaseInt Size() const { return cov_counters.size(); } - BaseInt operator[](ElementIndex i) const { - assert(i < ElementIndex(cov_counters.size())); + BaseInt operator[](IndexT i) const { + assert(i < IndexT(cov_counters.size())); return cov_counters[i]; } template BaseInt Cover(const IterableT& subset) { BaseInt covered = 0; - for (ElementIndex i : subset) { + for (IndexT i : subset) { covered += cov_counters[i] == 0 ? 1ULL : 0ULL; cov_counters[i]++; } @@ -56,7 +90,7 @@ class CoverCounters { template BaseInt Uncover(const IterableT& subset) { BaseInt uncovered = 0; - for (ElementIndex i : subset) { + for (IndexT i : subset) { --cov_counters[i]; uncovered += cov_counters[i] == 0 ? 1ULL : 0ULL; } @@ -67,19 +101,22 @@ class CoverCounters { template bool IsRedundantCover(IterableT const& subset) const { return absl::c_all_of(subset, - [&](ElementIndex i) { return cov_counters[i] > 0; }); + [&](IndexT i) { return cov_counters[i] > 0; }); } // Check if all the elements would still be covered if the subset was removed. template bool IsRedundantUncover(IterableT const& subset) const { return absl::c_all_of(subset, - [&](ElementIndex i) { return cov_counters[i] > 1; }); + [&](IndexT i) { return cov_counters[i] > 1; }); } private: - ElementToIntVector cov_counters; + util_intops::StrongVector cov_counters; }; + +using CoverCounters = CoverCountersImpl; +using FullCoverCounters = CoverCountersImpl; } // namespace Solution::Solution(const SubModel& model, @@ -106,9 +143,10 @@ BoundCBs::BoundCBs(const SubModel& model) prev_best_lb_(std::numeric_limits::lowest()), max_iter_countdown_(10 * model.num_focus_elements()), // Arbitrary from [1] - exit_test_countdown_(300), // Arbitrrary from [1] + exit_test_countdown_(300), // Arbitrary from [1] exit_test_period_(300), // Arbitrary from [1] - step_size_(0.1), // Arbitrary from [1] + unfixed_run_extension_(0), + step_size_(0.1), // Arbitrary from [1] last_min_lb_seen_(std::numeric_limits::max()), last_max_lb_seen_(.0), step_size_update_countdown_(20), // Arbitrary from [1] @@ -116,61 +154,71 @@ BoundCBs::BoundCBs(const SubModel& model) {} bool BoundCBs::ExitCondition(const SubgradientContext& context) { - if (last_core_update_countdown_ > 0) { - --last_core_update_countdown_; - return false; - } - if (--max_iter_countdown_ <= 0 || squared_norm_ <= kTol || - // Note: assumes integral costs - context.best_dual_state.lower_bound() >= - context.best_solution.cost() - context.model.fixed_cost() - .999) { + Cost best_lb = context.best_lower_bound; + Cost best_ub = context.best_solution.cost() - context.model.fixed_cost(); + if (--max_iter_countdown_ <= 0 || squared_norm_ <= kTol) { return true; } if (--exit_test_countdown_ > 0) { return false; } + if (prev_best_lb_ >= best_ub - CFT_BOUND_EPSILON) { + return true; + } exit_test_countdown_ = exit_test_period_; - Cost best_lower_bound = context.best_dual_state.lower_bound(); - Cost abs_improvement = best_lower_bound - prev_best_lb_; - Cost rel_improvement = DivideIfGE0(abs_improvement, best_lower_bound); - prev_best_lb_ = best_lower_bound; - return abs_improvement < 1.0 && rel_improvement < .001; + Cost abs_improvement = best_lb - prev_best_lb_; + Cost rel_improvement = DivideIfGE0(abs_improvement, best_lb); + prev_best_lb_ = best_lb; + + if (abs_improvement >= 1.0 || rel_improvement >= .001) { + return false; + } + + // (Not in [1]): During the first unfixed iteration we want to converge closer + // to the optimum + VLOG(3) << "[SUBG] First iteration extension " << unfixed_run_extension_; + BaseInt extension = (context.model.fixed_columns().empty() ? 2 : 1); + return unfixed_run_extension_++ >= extension; } void BoundCBs::ComputeMultipliersDelta(const SubgradientContext& context, ElementCostVector& delta_mults) { - Cost lower_bound = context.current_dual_state.lower_bound(); - UpdateStepSize(lower_bound); - direction_ = context.subgradient; - MakeMinimalCoverageSubgradient(context, direction_); - - if (prev_direction_.empty()) { - prev_direction_ = direction_; + BaseInt extension = (context.model.fixed_columns().empty() ? 2 : 1); + if (unfixed_run_extension_ < extension) { + MakeMinimalCoverageSubgradient(context, direction_); } squared_norm_ = .0; - // Stabilize current direction and compute squared norm for (ElementIndex i : context.model.ElementRange()) { - direction_[i] = stabilization_coeff * direction_[i] + - (1.0 - stabilization_coeff) * prev_direction_[i]; - if ((context.current_dual_state.multipliers()[i] <= .0) && - (direction_[i] < .0)) { + Cost curr_i_mult = context.current_dual_state.multipliers()[i]; + if ((curr_i_mult <= .0 && direction_[i] < .0) || + (curr_i_mult > CFT_MAX_MULTIPLIER && direction_[i] > .0)) { direction_[i] = 0; } + squared_norm_ += direction_[i] * direction_[i]; - prev_direction_[i] = direction_[i]; } + if (squared_norm_ <= kTol) { + delta_mults.assign(context.model.num_elements(), .0); + return; + } + + UpdateStepSize(context); Cost upper_bound = context.best_solution.cost() - context.model.fixed_cost(); + Cost lower_bound = context.current_dual_state.lower_bound(); Cost delta = upper_bound - lower_bound; - Cost step_constant = step_size_ * delta / squared_norm_; + Cost step_constant = (step_size_ * delta) / squared_norm_; + for (ElementIndex i : context.model.ElementRange()) { delta_mults[i] = step_constant * direction_[i]; + DCHECK(std::isfinite(delta_mults[i])); } } -void BoundCBs::UpdateStepSize(Cost lower_bound) { +void BoundCBs::UpdateStepSize(SubgradientContext context) { + Cost lower_bound = context.current_dual_state.lower_bound(); last_min_lb_seen_ = std::min(last_min_lb_seen_, lower_bound); last_max_lb_seen_ = std::max(last_max_lb_seen_, lower_bound); @@ -181,12 +229,10 @@ void BoundCBs::UpdateStepSize(Cost lower_bound) { Cost gap = DivideIfGE0(delta, last_max_lb_seen_); if (gap <= .001) { // Arbitray from [1] step_size_ *= 1.5; // Arbitray from [1] - - // Arbitrary from c4v4 - stabilization_coeff = (1.0 + stabilization_coeff) / 2.0; - + VLOG(4) << "[SUBG] Sep size set at " << step_size_; } else if (gap > .01) { // Arbitray from [1] step_size_ /= 2.0; // Arbitray from [1] + VLOG(4) << "[SUBG] Sep size set at " << step_size_; } last_min_lb_seen_ = std::numeric_limits::max(); last_max_lb_seen_ = .0; @@ -217,15 +263,15 @@ void BoundCBs::MakeMinimalCoverageSubgradient(const SubgradientContext& context, } } -bool BoundCBs::UpdateCoreModel(SubModel& core_model, - PrimalDualState& best_state) { - if (core_model.UpdateCore(best_state)) { - // grant at least 10 iterations before the next exit test - prev_best_lb_ = - std::min(prev_best_lb_, best_state.dual_state.lower_bound()); - exit_test_countdown_ = std::max(exit_test_countdown_, 20); - max_iter_countdown_ = std::max(max_iter_countdown_, 20); - last_core_update_countdown_ = 20; +bool BoundCBs::UpdateCoreModel(SubgradientContext context, + CoreModel& core_model, bool force) { + if (core_model.UpdateCore(context.best_lower_bound, context.best_multipliers, + context.best_solution, force)) { + prev_best_lb_ = std::numeric_limits::lowest(); + // Grant at least `min_iters` iterations before the next exit test + constexpr BaseInt min_iters = 10; + exit_test_countdown_ = std::max(exit_test_countdown_, min_iters); + max_iter_countdown_ = std::max(max_iter_countdown_, min_iters); return true; } return false; @@ -233,20 +279,33 @@ bool BoundCBs::UpdateCoreModel(SubModel& core_model, void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, PrimalDualState& best_state) { + CFT_MEASURE_SCOPE_DURATION(subgradient_time); DCHECK(ValidateSubModel(model)); ElementCostVector subgradient = ElementCostVector(model.num_elements(), .0); DualState dual_state = best_state.dual_state; + Cost best_lower_bound = dual_state.lower_bound(); + ElementCostVector best_multipliers = dual_state.multipliers(); Solution solution; ElementCostVector multipliers_delta(model.num_elements()); // to avoid allocs SubgradientContext context = {.model = model, .current_dual_state = dual_state, - .best_dual_state = best_state.dual_state, + .best_lower_bound = best_lower_bound, + .best_multipliers = best_multipliers, .best_solution = best_state.solution, .subgradient = subgradient}; - size_t iter = 0; - while (!cbs.ExitCondition(context)) { + + for (size_t iter = 1; !cbs.ExitCondition(context); ++iter) { + // Poor multipliers can lead to wasted iterations or stagnation in the + // subgradient method. To address this, we adjust the multipliers to + // get closer the trivial lower bound (= 0). + if (dual_state.lower_bound() < .0) { + VLOG(4) << "[SUBG] Dividing multipliers by 10"; + dual_state.DualUpdate( + model, [&](ElementIndex i, Cost& i_mult) { i_mult /= 10.0; }); + } + // Compute subgradient (O(nnz)) subgradient.assign(model.num_elements(), 1.0); for (SubsetIndex j : model.SubsetRange()) { @@ -259,10 +318,12 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, cbs.ComputeMultipliersDelta(context, multipliers_delta); dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { - i_mult = std::max(.0, i_mult + multipliers_delta[i]); + i_mult = std::clamp(i_mult + multipliers_delta[i], .0, + CFT_MAX_MULTIPLIER); }); - if (dual_state.lower_bound() > best_state.dual_state.lower_bound()) { - best_state.dual_state = dual_state; + if (dual_state.lower_bound() > best_lower_bound) { + best_lower_bound = dual_state.lower_bound(); + best_multipliers = dual_state.multipliers(); } cbs.RunHeuristic(context, solution); @@ -271,23 +332,35 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, best_state.solution = solution; } - if (++iter % 50 == 0) - std::cout << "Subgradient " << iter << " -- Bounds: Lower " - << dual_state.lower_bound() << ", best " - << best_state.dual_state.lower_bound() << " - Upper " - << best_state.solution.cost() - model.fixed_cost() - << ", global " << best_state.solution.cost() << "\n"; + VLOG_EVERY_N(4, 100) << "[SUBG] " << iter << ": Bounds: Lower " + << dual_state.lower_bound() << ", best " + << best_lower_bound << " - Upper " + << best_state.solution.cost() - model.fixed_cost() + << ", global " << best_state.solution.cost(); - if (cbs.UpdateCoreModel(model, best_state)) { - std::cout << "Core model has been updated.\n"; - dual_state = best_state.dual_state; + if (cbs.UpdateCoreModel(context, model)) { + dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { + i_mult = best_multipliers[i]; + }); + best_lower_bound = dual_state.lower_bound(); } } - std::cout << "Subgradient End" << " -- Bounds: Lower " - << dual_state.lower_bound() << ", best " - << best_state.dual_state.lower_bound() << " - Upper " - << best_state.solution.cost() - model.fixed_cost() << ", global " - << best_state.solution.cost() << "\n"; + + if (cbs.UpdateCoreModel(context, model, /*force=*/true)) { + dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { + i_mult = best_multipliers[i]; + }); + best_lower_bound = dual_state.lower_bound(); + } + best_state.dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { + i_mult = best_multipliers[i]; + }); + DCHECK_EQ(best_state.dual_state.lower_bound(), best_lower_bound); + + VLOG(3) << "[SUBG] End - Bounds: Lower " << dual_state.lower_bound() + << ", best " << best_lower_bound << " - Upper " + << best_state.solution.cost() - model.fixed_cost() << ", global " + << best_state.solution.cost(); } //////////////////////////////////////////////////////////////////////// @@ -544,6 +617,8 @@ Solution RunMultiplierBasedGreedy(const SubModel& model, Cost CoverGreedly(const SubModel& model, const DualState& dual_state, Cost cost_cutoff, BaseInt stop_size, std::vector& sol_subsets) { + CFT_MEASURE_SCOPE_DURATION(greedy_time); + Cost sol_cost = .0; for (SubsetIndex j : sol_subsets) { sol_cost += model.subset_costs()[j]; @@ -588,6 +663,7 @@ Cost CoverGreedly(const SubModel& model, const DualState& dual_state, // Either remove redundant columns or discard solution RedundancyRemover remover(model, covered_rows); // TODO(?): cache it! + return remover.TryRemoveRedundantCols(model, cost_cutoff, sol_subsets); } @@ -648,12 +724,12 @@ void FixBestColumns(SubModel& model, PrimalDualState& state) { fix_at_least, cols_to_fix); // Fix columns and update the model - Cost fixed_cost_delta = model.FixColumns(cols_to_fix); + Cost fixed_cost_delta = model.FixMoreColumns(cols_to_fix); - std::cout << "Fixed " << cols_to_fix.size() - << " new columns with cost: " << fixed_cost_delta << '\n'; - std::cout << "Globally fixed " << model.fixed_columns().size() - << " columns, with cost " << model.fixed_cost() << '\n'; + VLOG(3) << "[3FIX] Fixed " << cols_to_fix.size() + << " new columns with cost: " << fixed_cost_delta; + VLOG(3) << "[3FIX] Globally fixed " << model.fixed_columns().size() + << " columns, with cost " << model.fixed_cost(); // Update multipliers for the reduced model dual_state.DualUpdate(model, [&](ElementIndex core_i, Cost& multiplier) { @@ -665,6 +741,9 @@ void FixBestColumns(SubModel& model, PrimalDualState& state) { void RandomizeDualState(const SubModel& model, DualState& dual_state, std::mt19937& rnd) { + // In [1] this step is described, not completely sure if it actually helps + // or not. Seems to me one of those "throw in some randomness, it never hurts" + // thing. dual_state.DualUpdate(model, [&](ElementIndex, Cost& i_multiplier) { i_multiplier *= absl::Uniform(rnd, 0.9, 1.1); }); @@ -696,6 +775,7 @@ void HeuristicCBs::ComputeMultipliersDelta(const SubgradientContext& context, } PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { + CFT_MEASURE_SCOPE_DURATION(three_phase_time); DCHECK(ValidateSubModel(model)); PrimalDualState best_state = {.solution = init_solution, @@ -704,52 +784,188 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { best_state.solution = RunMultiplierBasedGreedy(model, best_state.dual_state); } - std::cout << "Initial lower bound: " << best_state.dual_state.lower_bound() - << "\nInitial solution cost: " << best_state.solution.cost() - << "\nStarting 3-phase algorithm\n"; + VLOG(2) << "[3PHS] Initial lower bound: " + << best_state.dual_state.lower_bound() + << ", Initial solution cost: " << best_state.solution.cost() + << ", Starting 3-phase algorithm"; PrimalDualState curr_state = best_state; BaseInt iter_count = 0; std::mt19937 rnd(0xcf7); - while (model.num_focus_elements() > 0 && - // note: assumes integral costs - curr_state.dual_state.lower_bound() < - best_state.solution.cost() - model.fixed_cost() - .999) { + while (model.num_focus_elements() > 0) { ++iter_count; - std::cout << "\n\n3Phase iteration: " << iter_count << '\n'; - std::cout << " Active size: rows " << model.num_focus_elements() << "/" - << model.num_elements() << ", columns " - << model.num_focus_subsets() << "/" << model.num_subsets() - << '\n'; + VLOG(2) << "[3PHS] 3Phase iteration: " << iter_count; + VLOG(2) << "[3PHS] Active size: rows " << model.num_focus_elements() << "/" + << model.num_elements() << ", columns " << model.num_focus_subsets() + << "/" << model.num_subsets(); // Phase 1: refine the current dual_state and model BoundCBs dual_bound_cbs(model); - std::cout << "\nSubgradient Phase:\n"; + VLOG(2) << "[3PHS] Subgradient Phase:"; SubgradientOptimization(model, dual_bound_cbs, curr_state); if (iter_count == 1) { best_state.dual_state = curr_state.dual_state; } + if (curr_state.dual_state.lower_bound() >= + best_state.solution.cost() - model.fixed_cost() - CFT_BOUND_EPSILON) { + break; + } + // Phase 2: search for good solutions HeuristicCBs heuristic_cbs; heuristic_cbs.set_step_size(dual_bound_cbs.step_size()); - std::cout << "\nHeuristic Phase:\n"; + VLOG(2) << "[3PHS] Heuristic Phase:"; SubgradientOptimization(model, heuristic_cbs, curr_state); + if (iter_count == 1 && best_state.dual_state.lower_bound() < + curr_state.dual_state.lower_bound()) { + best_state.dual_state = curr_state.dual_state; + } if (curr_state.solution.cost() < best_state.solution.cost()) { best_state.solution = curr_state.solution; } + if (curr_state.dual_state.lower_bound() >= + best_state.solution.cost() - model.fixed_cost() - CFT_BOUND_EPSILON) { + break; + } - std::cout << "\n3Phase Bounds: Lower " - << best_state.dual_state.lower_bound() << ", Upper " - << best_state.solution.cost() << '\n'; + VLOG(2) << "[3PHS] 3Phase Bounds: Residual Lower " + << curr_state.dual_state.lower_bound() + model.fixed_cost() + << ", Upper " << best_state.solution.cost(); // Phase 3: Fix the best columns (diving) FixBestColumns(model, curr_state); RandomizeDualState(model, curr_state.dual_state, rnd); } - std::cout << "\n\n3Phase End\n"; - std::cout << "Final Bounds: Lower " << best_state.dual_state.lower_bound() - << ", Upper " << best_state.solution.cost() << '\n'; + VLOG(2) << "[3PHS] 3Phase End: "; + VLOG(2) << "[3PHS] Bounds: Residual Lower " + << curr_state.dual_state.lower_bound() + model.fixed_cost() + << ", Upper " << best_state.solution.cost(); + + return best_state; +} + +/////////////////////////////////////////////////////////////////////// +///////////////////// OUTER REFINEMENT PROCEDURE ////////////////////// +/////////////////////////////////////////////////////////////////////// + +namespace { + +struct GapContribution { + Cost gap; + FullSubsetIndex idx; +}; + +std::vector SelectColumnByGapContribution( + const SubModel& model, const PrimalDualState& best_state, + BaseInt nrows_to_fix) { + const auto& [solution, dual_state] = best_state; + + FullCoverCounters row_coverage(model.num_elements()); + auto full_model = model.StrongTypedFullModelView(); + + for (FullSubsetIndex j : solution.subsets()) { + row_coverage.Cover(full_model.columns()[j]); + } + + std::vector gap_contributions; + for (FullSubsetIndex j : solution.subsets()) { + Cost j_gap = .0; + Cost reduced_cost = dual_state.reduced_costs()[static_cast(j)]; + for (FullElementIndex i : full_model.columns()[j]) { + Cost i_mult = dual_state.multipliers()[static_cast(i)]; + j_gap += i_mult * (1.0 - 1.0 / row_coverage[i]); + reduced_cost -= i_mult; + } + j_gap += std::max(reduced_cost, 0.0); + gap_contributions.push_back({j_gap, j}); + } + absl::c_sort(gap_contributions, [](GapContribution g1, GapContribution g2) { + return g1.gap < g2.gap; + }); + + BaseInt covered_rows = 0; + row_coverage.Reset(model.num_elements()); + std::vector cols_to_fix; + for (auto [j_gap, j] : gap_contributions) { + covered_rows += row_coverage.Cover(full_model.columns()[j]); + if (covered_rows > nrows_to_fix) { + break; + } + cols_to_fix.push_back(static_cast(j)); + } + return cols_to_fix; +} +} // namespace + +PrimalDualState RunCftHeuristic(SubModel& model, + const Solution& init_solution) { + CFT_MEASURE_SCOPE_DURATION(refinement_time); + + PrimalDualState best_state = {.solution = init_solution, + .dual_state = MakeTentativeDualState(model)}; + if (best_state.solution.Empty()) { + best_state.solution = + RunMultiplierBasedGreedy(model, best_state.dual_state); + } + + Cost cost_cutoff = std::numeric_limits::lowest(); + double fix_minimum = .3; // Arbitrary from [1] + double fix_increment = 1.1; // Arbitrary from [1] + double fix_fraction = fix_minimum; + + for (BaseInt iter_counter = 0; model.num_elements() > 0; ++iter_counter) { + VLOG(1) << "[CFTH] Refinement iteration: " << iter_counter; + Cost fixed_cost_before = model.fixed_cost(); + auto [solution, dual_state] = RunThreePhase(model, best_state.solution); + if (iter_counter == 0) { + best_state.dual_state = std::move(dual_state); + } + if (solution.cost() < best_state.solution.cost()) { + best_state.solution = solution; + fix_fraction = fix_minimum; + } + cost_cutoff = std::max(cost_cutoff, + fixed_cost_before + dual_state.lower_bound()); + VLOG(1) << "[CFTH] Refinement Bounds: Residual Lower " << cost_cutoff + << ", Real Lower " << best_state.dual_state.lower_bound() + << ", Upper " << best_state.solution.cost(); + if (best_state.solution.cost() - CFT_BOUND_EPSILON <= cost_cutoff) { + break; + } + + fix_fraction = std::min(1.0, fix_fraction * fix_increment); + std::vector cols_to_fix = SelectColumnByGapContribution( + model, best_state, model.num_elements() * fix_fraction); + + if (!cols_to_fix.empty()) { + model.ResetColumnFixing(cols_to_fix, best_state.dual_state); + } + VLOG(1) << "[CFTH] Fixed " << cols_to_fix.size() + << " new columns with cost: " << model.fixed_cost(); + VLOG(1) << "[CFTH] Model sizes: rows " << model.num_focus_elements() << "/" + << model.num_elements() << ", columns " << model.num_focus_subsets() + << "/" << model.num_subsets(); + } + +#ifdef CFT_MEASURE_TIME + double subg_t = subgradient_time.Get(); + double greedy_t = greedy_time.Get(); + double three_phase_t = three_phase_time.Get(); + double refinement_t = refinement_time.Get(); + + printf("Subgradient time: %8.2f (%.1f%%)\n", subg_t, + 100 * subg_t / refinement_t); + printf("Greedy Heur time: %8.2f (%.1f%%)\n", greedy_t, + 100 * greedy_t / refinement_t); + printf("SubG - Greedy time: %8.2f (%.1f%%)\n", subg_t - greedy_t, + 100 * (subg_t - greedy_t) / refinement_t); + printf("3Phase time: %8.2f (%.1f%%)\n", three_phase_t, + 100 * three_phase_t / refinement_t); + printf("3Phase - Subg time: %8.2f (%.1f%%)\n", three_phase_t - subg_t, + 100 * (three_phase_t - subg_t) / refinement_t); + printf("Total CFT time: %8.2f (%.1f%%)\n", refinement_t, 100.0); +#endif return best_state; } @@ -760,9 +976,16 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { namespace { std::vector ComputeTentativeFocus(StrongModelView full_model) { - FullSubsetBoolVector selected(full_model.num_subsets(), false); std::vector columns_focus; + + if (full_model.num_subsets() <= 2 * kMinCov * full_model.num_elements()) { + columns_focus.resize(full_model.num_subsets()); + absl::c_iota(columns_focus, FullSubsetIndex(0)); + return columns_focus; + } + columns_focus.reserve(full_model.num_elements() * kMinCov); + FullSubsetBoolVector selected(full_model.num_subsets(), false); // Select the first min_row_coverage columns for each row for (const auto& row : full_model.rows()) { @@ -781,82 +1004,94 @@ std::vector ComputeTentativeFocus(StrongModelView full_model) { absl::c_sort(columns_focus); return columns_focus; } - -void SelecteMinRedCostColumns(FilterModelView full_model, - const SubsetCostVector& reduced_costs, - std::vector& new_core_columns, - FullSubsetBoolVector& selected) { - DCHECK_EQ(reduced_costs.size(), full_model.num_subsets()); - DCHECK_EQ(selected.size(), full_model.num_subsets()); - - std::vector candidates; - for (SubsetIndex j : full_model.SubsetRange()) - if (reduced_costs[j] < 0.1) { - candidates.push_back(j); - } - - BaseInt max_size = 5 * full_model.num_focus_elements(); - if (candidates.size() > max_size) { - absl::c_nth_element(candidates, candidates.begin() + max_size - 1, - [&](SubsetIndex j1, SubsetIndex j2) { - return reduced_costs[j1] < reduced_costs[j2]; - }); - candidates.resize(max_size); - } - for (SubsetIndex j : candidates) { - FullSubsetIndex j_full = static_cast(j); - if (!selected[j_full]) { - selected[j_full] = true; - new_core_columns.push_back(j_full); - } - } -} - -static void SelectMinRedCostByRow(FilterModelView full_model, - const SubsetCostVector& reduced_costs, - std::vector& columns_map, - FullSubsetBoolVector& selected) { - DCHECK_EQ(reduced_costs.size(), full_model.num_subsets()); - DCHECK_EQ(selected.size(), full_model.num_subsets()); - - for (const auto& row : full_model.rows()) { - // Collect best `kMinCov` columns covering row `i` - SubsetIndex best_cols[kMinCov]; - BaseInt best_size = 0; - for (SubsetIndex j : row) { - if (best_size < kMinCov) { - best_cols[best_size++] = j; - continue; - } - if (reduced_costs[j] < reduced_costs[best_cols[kMinCov - 1]]) { - BaseInt n = kMinCov - 1; - while (n > 0 && reduced_costs[j] < reduced_costs[best_cols[n - 1]]) { - best_cols[n] = best_cols[n - 1]; - --n; - } - best_cols[n] = j; - } - } - - DCHECK(best_size > 0); - for (BaseInt s = 0; s < best_size; ++s) { - FullSubsetIndex j = static_cast(best_cols[s]); - if (!selected[j]) { - selected[j] = true; - columns_map.push_back(j); - } - } - } -} } // namespace +const std::vector& +FullToCoreModel::ColumnSelector::ComputeNewSelection( + FilterModelView full_model, + const std::vector& forced_columns, + const SubsetCostVector& reduced_costs) { + selected_.assign(full_model.num_subsets(), false); + row_cover_counts_.assign(full_model.num_elements(), 0); + rows_left_to_cover_ = full_model.num_focus_elements(); + selection_.clear(); + candidates_.clear(); + + // Always retain best solution in the core model (if possible) + for (FullSubsetIndex full_j : forced_columns) { + SubsetIndex j = static_cast(full_j); + if (full_model.IsFocusCol(j)) { + SelectColumn(full_model, j); + } + } + + auto subset_range = full_model.SubsetRange(); + candidates_ = {subset_range.begin(), subset_range.end()}; + absl::c_sort(candidates_, [&](auto j1, auto j2) { + return reduced_costs[j1] < reduced_costs[j2]; + }); + first_unselected_ = candidates_.begin(); + + SelecteMinRedCostColumns(full_model, reduced_costs); + SelectMinRedCostByRow(full_model, reduced_costs); + absl::c_sort(selection_); + return selection_; +} + +bool FullToCoreModel::ColumnSelector::SelectColumn(FilterModelView full_model, + SubsetIndex j) { + if (selected_[j]) { + return false; + } + for (ElementIndex i : full_model.columns()[j]) { + selected_[j] = true; // Detect empty columns + if (++row_cover_counts_[i] == kMinCov) { + --rows_left_to_cover_; + } + } + if (selected_[j]) { // Skip empty comlumns + selection_.push_back(static_cast(j)); + } + return selected_[j]; +} + +void FullToCoreModel::ColumnSelector::SelecteMinRedCostColumns( + FilterModelView full_model, const SubsetCostVector& reduced_costs) { + BaseInt selected_size = 0; + BaseInt max_size = kMinCov * full_model.num_elements(); + auto it = first_unselected_; + while (it != candidates_.end() && reduced_costs[*it] < 0.1 && + selected_size < max_size) { + selected_size += SelectColumn(full_model, *it++) ? 1 : 0; + } + first_unselected_ = it; +} + +void FullToCoreModel::ColumnSelector::SelectMinRedCostByRow( + FilterModelView full_model, const SubsetCostVector& reduced_costs) { + auto it = first_unselected_; + while (it != candidates_.end() && rows_left_to_cover_ > 0) { + SubsetIndex j = *it++; + if (selected_[j]) { + continue; + } + for (ElementIndex i : full_model.columns()[j]) { + ++row_cover_counts_[i]; + rows_left_to_cover_ += (row_cover_counts_[i] == kMinCov ? 1 : 0); + selected_[j] = selected_[j] || (row_cover_counts_[i] <= kMinCov); + } + if (selected_[j]) { + selection_.push_back(static_cast(j)); + } + } +} + FullToCoreModel::FullToCoreModel(const Model* full_model) : SubModel(full_model, ComputeTentativeFocus(StrongModelView(full_model))), full_model_(full_model), is_focus_col_(full_model->num_subsets(), true), is_focus_row_(full_model->num_elements(), true), - num_subsets_(full_model->num_subsets()), - num_elements_(full_model->num_elements()), + prev_best_lower_bound_(std::numeric_limits::lowest()), full_dual_state_(*full_model), best_dual_state_(full_dual_state_) { ResetPricingPeriod(); @@ -870,15 +1105,19 @@ void FullToCoreModel::ResetPricingPeriod() { update_max_period_ = std::min(1000, full_model_->num_elements() / 3); } -Cost FullToCoreModel::FixColumns( +Cost FullToCoreModel::FixMoreColumns( const std::vector& columns_to_fix) { StrongModelView typed_full_model = StrongTypedFullModelView(); for (SubsetIndex core_j : columns_to_fix) { FullSubsetIndex full_j = SubModel::MapCoreToFullSubsetIndex(core_j); + DCHECK(IsFocusCol(full_j)); IsFocusCol(full_j) = false; + bool any_active = false; for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { + any_active |= IsFocusRow(full_i); IsFocusRow(full_i) = false; } + DCHECK(any_active); } for (FullSubsetIndex full_j : typed_full_model.SubsetRange()) { if (!IsFocusCol(full_j)) { @@ -892,59 +1131,95 @@ Cost FullToCoreModel::FixColumns( } } } + ResetPricingPeriod(); - Cost fixed_cost = base::FixColumns(columns_to_fix); + Cost fixed_cost = base::FixMoreColumns(columns_to_fix); DCHECK(FullToSubModelInvariantCheck()); return fixed_cost; } -bool FullToCoreModel::UpdateCore(PrimalDualState& core_state) { - if (--update_countdown_ > 0) { - return false; - } +void FullToCoreModel::ResetColumnFixing( + const std::vector& full_columns_to_fix, + const DualState& dual_state) { + is_focus_col_.assign(full_model_->num_subsets(), true); + is_focus_row_.assign(full_model_->num_elements(), true); - UpdateDualState(core_state.dual_state, full_dual_state_, best_dual_state_); - UpdatePricingPeriod(full_dual_state_, core_state); - std::cout << "Lower bounds: Real " << full_dual_state_.lower_bound() - << ", Core " << core_state.dual_state.lower_bound() << '\n'; + full_dual_state_ = dual_state; - auto fixing_full_model = FixingFullModelView(); - FullSubsetBoolVector selected_columns(fixing_full_model.num_subsets(), false); - std::vector new_core_columns; + // We could implement and in-place core-model update that removes old fixings, + // set the new one while also updating the column focus. This solution is much + // simpler. It just create a new core-model object from scratch and then uses + // the existing interface. + const auto& selection = col_selector_.ComputeNewSelection( + FixingFullModelView(), full_columns_to_fix, + full_dual_state_.reduced_costs()); + static_cast(*this) = SubModel(full_model_, selection); - // Always retain best solution in the core model - for (FullSubsetIndex full_j : core_state.solution.subsets()) { - if (IsFocusCol(full_j)) { - new_core_columns.push_back(full_j); - selected_columns[full_j] = true; + // TODO(anyone): Improve this. It's Inefficient but hardly a botleneck and it + // also avoid storing a full->core column map. + std::vector columns_to_fix; + for (SubsetIndex core_j : SubsetRange()) { + for (FullSubsetIndex full_j : full_columns_to_fix) { + if (full_j == MapCoreToFullSubsetIndex(core_j)) { + columns_to_fix.push_back(core_j); + break; + } } } + DCHECK_EQ(columns_to_fix.size(), full_columns_to_fix.size()); + FixMoreColumns(columns_to_fix); + DCHECK(FullToSubModelInvariantCheck()); +} - SelecteMinRedCostColumns(fixing_full_model, full_dual_state_.reduced_costs(), - new_core_columns, selected_columns); - SelectMinRedCostByRow(fixing_full_model, full_dual_state_.reduced_costs(), - new_core_columns, selected_columns); +void FullToCoreModel::SizeUpdate() { + is_focus_col_.resize(full_model_->num_subsets(), true); +} - // NOTE: unnecessary, but it keeps equivalence between SubModelView/SubModel - absl::c_sort(new_core_columns); - SetFocus(new_core_columns); - core_state.dual_state.DualUpdate(*this, [](ElementIndex i, Cost& i_mult) { - // multipliers didn't cange, but reduced cost must be recomputed - }); +bool FullToCoreModel::IsTimeToUpdate(Cost best_lower_bound, bool force) { + if (!force && --update_countdown_ > 0) { + return false; + } + if (best_lower_bound == prev_best_lower_bound_) { + return false; + } + prev_best_lower_bound_ = best_lower_bound; + return true; +} +void FullToCoreModel::ComputeAndSetFocus(Cost best_lower_bound, + const Solution& best_solution) { + const auto& selection = col_selector_.ComputeNewSelection( + FixingFullModelView(), best_solution.subsets(), + full_dual_state_.reduced_costs()); + base::SetFocus(selection); + UpdatePricingPeriod(full_dual_state_, best_lower_bound, + best_solution.cost() - fixed_cost()); + VLOG(3) << "[F2CU] Core-update: Lower bounds: real " + << full_dual_state_.lower_bound() << ", core " << best_lower_bound + << ", core size: " << num_focus_elements() << "x" + << num_focus_subsets(); +} + +bool FullToCoreModel::UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) { + if (!IsTimeToUpdate(best_lower_bound, force)) { + return false; + } + UpdateMultipliers(best_multipliers); + ComputeAndSetFocus(best_lower_bound, best_solution); DCHECK(FullToSubModelInvariantCheck()); return true; } void FullToCoreModel::UpdatePricingPeriod(const DualState& full_dual_state, - const PrimalDualState& core_state) { - DCHECK_GE(core_state.dual_state.lower_bound() + 1e-6, - full_dual_state.lower_bound()); - DCHECK_GE(core_state.solution.cost(), .0); + Cost core_lower_bound, + Cost core_upper_bound) { + DCHECK_GE(core_lower_bound + 1e-6, full_dual_state.lower_bound()); + DCHECK_GE(core_upper_bound, .0); - Cost delta = - core_state.dual_state.lower_bound() - full_dual_state.lower_bound(); - Cost ratio = DivideIfGE0(delta, core_state.solution.cost()); + Cost delta = core_lower_bound - full_dual_state.lower_bound(); + Cost ratio = DivideIfGE0(delta, core_upper_bound); if (ratio <= 1e-6) { update_period_ = std::min(update_max_period_, 10 * update_period_); } else if (ratio <= 0.02) { @@ -957,15 +1232,14 @@ void FullToCoreModel::UpdatePricingPeriod(const DualState& full_dual_state, update_countdown_ = update_period_; } -void FullToCoreModel::UpdateDualState(const DualState& core_dual_state, - DualState& full_dual_state, - DualState& best_dual_state) { +Cost FullToCoreModel::UpdateMultipliers( + const ElementCostVector& core_multipliers) { auto fixing_full_model = FixingFullModelView(); full_dual_state_.DualUpdate( fixing_full_model, [&](ElementIndex full_i, Cost& i_mult) { ElementIndex core_i = MapFullToCoreElementIndex(static_cast(full_i)); - i_mult = core_dual_state.multipliers()[core_i]; + i_mult = core_multipliers[core_i]; }); // Here, we simply check if any columns have been fixed, and only update the @@ -983,6 +1257,7 @@ void FullToCoreModel::UpdateDualState(const DualState& core_dual_state, full_dual_state_.lower_bound() > best_dual_state_.lower_bound()) { best_dual_state_ = full_dual_state_; } + return full_dual_state_.lower_bound(); } bool FullToCoreModel::FullToSubModelInvariantCheck() { @@ -1009,8 +1284,50 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { " not found in full model view.\n"); return false; } + + const auto& core_column = sub_model.columns()[core_j]; + if (core_column.begin() == core_column.end()) { + std::cerr << absl::StrCat("Core subset ", core_j, " empty.\n"); + return false; + } + + const auto& full_column = typed_full_model.columns()[full_j]; + if (full_column.begin() == full_column.end()) { + std::cerr << absl::StrCat("Full subset ", full_j, " empty.\n"); + return false; + } + + // Assumes corresponding elements have the same order in both models + auto core_it = core_column.begin(); + for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { + if (sub_model.MapFullToCoreElementIndex(full_i) != *core_it) { + continue; + } + if (sub_model.MapCoreToFullElementIndex(*core_it) != full_i) { + std::cerr << absl::StrCat( + "Subset ", core_j, " in sub-model has mapped element ", *core_it, + " but it is not the same as the full model.\n"); + return false; + } + if (++core_it == core_column.end()) { + break; + } + } + if (core_it != core_column.end()) { + std::cerr << absl::StrCat("Subset ", core_j, + " in sub-model has no mapped element ", + *core_it, " in full model view.\n"); + return false; + } } + for (ElementIndex core_i : sub_model.ElementRange()) { + const auto& core_row = sub_model.rows()[core_i]; + if (core_row.begin() == core_row.end()) { + std::cerr << absl::StrCat("Core row ", core_i, " empty.\n"); + return false; + } + FullElementIndex full_i = sub_model.MapCoreToFullElementIndex(core_i); if (!is_focus_row_[static_cast(full_i)]) { std::cerr << absl::StrCat("Element ", core_i, @@ -1019,6 +1336,7 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { return false; } } + for (FullElementIndex full_i : typed_full_model.ElementRange()) { if (!IsFocusRow(full_i)) { continue; diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 5cdc99123e..9481161a57 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -100,7 +100,7 @@ class Solution { } private: - Cost cost_; + Cost cost_ = std::numeric_limits::max(); std::vector subsets_; }; @@ -122,8 +122,9 @@ inline Cost DivideIfGE0(Cost numerator, Cost denominator) { class DualState { public: DualState() = default; + DualState(const DualState&) = default; template - DualState(const SubModelT& model) + explicit DualState(const SubModelT& model) : lower_bound_(.0), multipliers_(model.num_elements(), .0), reduced_costs_(model.subset_costs().begin(), @@ -143,6 +144,7 @@ class DualState { for (ElementIndex i : model.ElementRange()) { multiplier_operator(i, multipliers_[i]); lower_bound_ += multipliers_[i]; + DCHECK(std::isfinite(multipliers_[i])); DCHECK_GE(multipliers_[i], .0); } lower_bound_ += ComputeReducedCosts(model, multipliers_, reduced_costs_); @@ -188,7 +190,11 @@ struct PrimalDualState { struct SubgradientContext { const SubModel& model; const DualState& current_dual_state; - const DualState& best_dual_state; + + // Avoid copying unused reduced cost during subgradient + const Cost& best_lower_bound; + const ElementCostVector& best_multipliers; + const Solution& best_solution; const ElementCostVector& subgradient; }; @@ -201,7 +207,8 @@ class SubgradientCBs { virtual void RunHeuristic(const SubgradientContext&, Solution&) = 0; virtual void ComputeMultipliersDelta(const SubgradientContext&, ElementCostVector& delta_mults) = 0; - virtual bool UpdateCoreModel(SubModel&, PrimalDualState&) = 0; + virtual bool UpdateCoreModel(SubgradientContext context, + CoreModel& core_model, bool force = false) = 0; virtual ~SubgradientCBs() = default; }; @@ -223,8 +230,8 @@ class BoundCBs : public SubgradientCBs { ElementCostVector& delta_mults) override; void RunHeuristic(const SubgradientContext& context, Solution& solution) override {} - bool UpdateCoreModel(SubModel& core_model, - PrimalDualState& best_state) override; + bool UpdateCoreModel(SubgradientContext context, CoreModel& core_model, + bool force = false) override; private: void MakeMinimalCoverageSubgradient(const SubgradientContext& context, @@ -232,24 +239,7 @@ class BoundCBs : public SubgradientCBs { private: Cost squared_norm_; - - // This addition implements a simplified stabilization technique inspired by - // works in the literature, such as: - // Frangioni, A., Gendron, B., Gorgone, E. (2015): - // "On the computational efficiency of subgradient methods: A case study in - // combinatorial optimization." - // - // The approach aims to reduce oscillations (zig-zagging) in the subgradient - // by using a "moving average" of the current and previous subgradients. The - // current subgradient is weighted by a factor of alpha, while the previous - // subgradients contribution is weighted by (1 - alpha). The parameter alpha - // is set to 0.5 by default but can be adjusted for tuning. The resulting - // stabilized subgradient is referred to as "direction" to distinguish it from - // the original subgradient. - Cost stabilization_coeff = 0.5; // Arbitrary from c4v4 ElementCostVector direction_; - ElementCostVector prev_direction_; - std::vector lagrangian_solution_; // Stopping condition @@ -257,10 +247,10 @@ class BoundCBs : public SubgradientCBs { BaseInt max_iter_countdown_; BaseInt exit_test_countdown_; BaseInt exit_test_period_; - BaseInt last_core_update_countdown_; + BaseInt unfixed_run_extension_; // Step size - void UpdateStepSize(Cost lower_bound); + void UpdateStepSize(SubgradientContext context); Cost step_size_; Cost last_min_lb_seen_; Cost last_max_lb_seen_; @@ -301,14 +291,15 @@ class HeuristicCBs : public SubgradientCBs { bool ExitCondition(const SubgradientContext& context) override { Cost upper_bound = context.best_solution.cost() - context.model.fixed_cost(); - Cost lower_bound = context.best_dual_state.lower_bound(); + Cost lower_bound = context.best_lower_bound; return upper_bound - .999 < lower_bound || --countdown_ <= 0; } void RunHeuristic(const SubgradientContext& context, Solution& solution) override; void ComputeMultipliersDelta(const SubgradientContext& context, ElementCostVector& delta_mults) override; - bool UpdateCoreModel(SubModel& model, PrimalDualState& state) override { + bool UpdateCoreModel(SubgradientContext context, CoreModel& core_model, + bool force = false) override { return false; } @@ -320,10 +311,20 @@ class HeuristicCBs : public SubgradientCBs { PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution = {}); +/////////////////////////////////////////////////////////////////////// +///////////////////// OUTER REFINEMENT PROCEDURE ////////////////////// +/////////////////////////////////////////////////////////////////////// + +PrimalDualState RunCftHeuristic(SubModel& model, + const Solution& init_solution = {}); + /////////////////////////////////////////////////////////////////////// //////////////////////// FULL TO CORE PRICING ///////////////////////// /////////////////////////////////////////////////////////////////////// +// Coverage counter to decide the number of columns to keep in the core model. +static constexpr BaseInt kMinCov = 5; + // CoreModel extractor. Stores a pointer to the full model and specilized // UpdateCore in such a way to updated the SubModel (stored as base class) and // focus the search on a small windows of the full model. @@ -335,14 +336,63 @@ class FullToCoreModel : public SubModel { BaseInt max_period; }; + // This class handles the logic for selecting columns based on their reduced + // costs and the number of rows they cover. While this implementation is more + // complex than what would typically be required for static `SetCoverModel`s, + // it is designed to efficiently handle dynamically updated models where new + // columns are generated over time. In this case, recomputing the row view + // from scratch each time would introduce significant overhead. To avoid this, + // the column selection logic operates solely on the column view, without + // relying on the row view. + // + // NOTE: A cleaner alternative would involve modifying the `SetCoverModel` + // implementation to support incremental updates to the row view as new + // columns are added. This approach would reduce overhead while enabling a + // simpler and more efficient column selection process. + // + // NOTE: The row-view based approach is available at commit: + // a598cf83930629853f72b964ebcff01f7a9378e0 + class ColumnSelector { + public: + const std::vector& ComputeNewSelection( + FilterModelView full_model, + const std::vector& forced_columns, + const SubsetCostVector& reduced_costs); + + private: + bool SelectColumn(FilterModelView full_model, SubsetIndex j); + void SelecteMinRedCostColumns(FilterModelView full_model, + const SubsetCostVector& reduced_costs); + void SelectMinRedCostByRow(FilterModelView full_model, + const SubsetCostVector& reduced_costs); + + private: + std::vector candidates_; + std::vector::const_iterator first_unselected_; + ElementToIntVector row_cover_counts_; + BaseInt rows_left_to_cover_; + + std::vector selection_; + SubsetBoolVector selected_; + }; + public: + FullToCoreModel() = default; FullToCoreModel(const Model* full_model); - Cost FixColumns(const std::vector& columns_to_fix) override; - bool UpdateCore(PrimalDualState& core_state) override; + Cost FixMoreColumns(const std::vector& columns_to_fix) override; + void ResetColumnFixing(const std::vector& columns_to_fix, + const DualState& state) override; + bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) override; void ResetPricingPeriod(); const DualState& best_dual_state() const { return best_dual_state_; } + bool FullToSubModelInvariantCheck(); + + protected: + void SizeUpdate(); + bool IsTimeToUpdate(Cost best_lower_bound, bool force); - private: decltype(auto) IsFocusCol(FullSubsetIndex j) { return is_focus_col_[static_cast(j)]; } @@ -350,9 +400,9 @@ class FullToCoreModel : public SubModel { return is_focus_row_[static_cast(i)]; } void UpdatePricingPeriod(const DualState& full_dual_state, - const PrimalDualState& core_state); - void UpdateDualState(const DualState& core_dual_state, - DualState& full_dual_state, DualState& best_dual_state); + Cost core_lower_bound, Cost core_upper_bound); + Cost UpdateMultipliers(const ElementCostVector& core_multipliers); + void ComputeAndSetFocus(Cost best_lower_bound, const Solution& best_solution); // Views are not composable (for now), so we can either access the full_model // with the strongly typed view or with the filtered view. @@ -360,7 +410,8 @@ class FullToCoreModel : public SubModel { // Access the full model filtered by the current columns fixed. FilterModelView FixingFullModelView() const { return FilterModelView(full_model_, &is_focus_col_, &is_focus_row_, - num_subsets_, num_elements_); + full_model_->num_subsets(), + full_model_->num_elements()); } // Access the full model with the strongly typed view. @@ -368,7 +419,8 @@ class FullToCoreModel : public SubModel { return StrongModelView(full_model_); } - bool FullToSubModelInvariantCheck(); + std::vector SelectNewCoreColumns( + const std::vector& forced_columns = {}); private: const Model* full_model_; @@ -388,19 +440,17 @@ class FullToCoreModel : public SubModel { // implementation that avoids this memory optimization was preferred. ElementBoolVector is_focus_row_; - BaseInt num_subsets_; - BaseInt num_elements_; - + BaseInt selection_coefficient_ = kMinCov; + Cost prev_best_lower_bound_; DualState full_dual_state_; DualState best_dual_state_; BaseInt update_countdown_; BaseInt update_period_; BaseInt update_max_period_; -}; -// Coverage counter to decide the number of columns to keep in the core model. -static constexpr BaseInt kMinCov = 5; + ColumnSelector col_selector_; // Here to avoid reallocations +}; } // namespace operations_research::scp diff --git a/ortools/set_cover/set_cover_submodel.cc b/ortools/set_cover/set_cover_submodel.cc index 02a13eb255..59202e9f92 100644 --- a/ortools/set_cover/set_cover_submodel.cc +++ b/ortools/set_cover/set_cover_submodel.cc @@ -14,37 +14,15 @@ #include "ortools/set_cover/set_cover_submodel.h" #include "ortools/base/stl_util.h" +#include "ortools/set_cover/set_cover_cft.h" +#include "ortools/set_cover/set_cover_views.h" namespace operations_research::scp { -template -void PrintSubModelSummary(const SubModelT& model) { - std::cout << "SubModelView sizes: rows " << model.num_focus_elements() << "/" - << model.num_elements() << ", columns " << model.num_focus_subsets() - << "/" << model.num_subsets() << '\n'; - std::cout << "Fixing: columns " << model.fixed_columns().size() << " cost " - << model.fixed_cost() << '\n'; -} - SubModelView::SubModelView(const Model* model) : base_view(model, &cols_sizes_, &rows_sizes_, &cols_focus_, &rows_focus_), - full_model_(model), - cols_sizes_(model->num_subsets()), - rows_sizes_(model->num_elements()) { - DCHECK(full_model_ != nullptr); - cols_focus_.reserve(model->num_subsets()); - rows_focus_.reserve(model->num_elements()); - for (SubsetIndex j : model->SubsetRange()) { - cols_sizes_[j] = model->columns()[j].size(); - cols_focus_.push_back(j); - } - for (ElementIndex i : model->ElementRange()) { - rows_sizes_[i] = model->rows()[i].size(); - rows_focus_.push_back(i); - } - fixed_columns_.clear(); - fixed_cost_ = 0.0; - PrintSubModelSummary(*this); + full_model_(model) { + ResetToIdentitySubModel(); DCHECK(ValidateSubModel(*this)); } @@ -59,7 +37,25 @@ SubModelView::SubModelView(const Model* model, SetFocus(columns_focus); } -Cost SubModelView::FixColumns(const std::vector& columns_to_fix) { +void SubModelView::ResetToIdentitySubModel() { + cols_sizes_.resize(full_model_->num_subsets()); + rows_sizes_.resize(full_model_->num_elements()); + cols_focus_.clear(); + rows_focus_.clear(); + for (SubsetIndex j : full_model_->SubsetRange()) { + cols_sizes_[j] = full_model_->columns()[j].size(); + cols_focus_.push_back(j); + } + for (ElementIndex i : full_model_->ElementRange()) { + rows_sizes_[i] = full_model_->rows()[i].size(); + rows_focus_.push_back(i); + } + fixed_columns_.clear(); + fixed_cost_ = .0; +} + +Cost SubModelView::FixMoreColumns( + const std::vector& columns_to_fix) { DCHECK(full_model_ != nullptr); Cost old_fixed_cost = fixed_cost_; if (columns_to_fix.empty()) { @@ -88,11 +84,20 @@ Cost SubModelView::FixColumns(const std::vector& columns_to_fix) { gtl::STLEraseAllFromSequenceIf( &rows_focus_, [&](ElementIndex i) { return !rows_sizes_[i]; }); - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); return fixed_cost_ - old_fixed_cost; } +void SubModelView::ResetColumnFixing( + const std::vector& columns_to_fix, const DualState&) { + ResetToIdentitySubModel(); + std::vector core_column_to_fix; + for (FullSubsetIndex full_j : columns_to_fix) { + core_column_to_fix.push_back(static_cast(full_j)); + } + FixMoreColumns(core_column_to_fix); +} + void SubModelView::SetFocus(const std::vector& columns_focus) { DCHECK(full_model_ != nullptr); DCHECK(!rows_sizes_.empty()); @@ -125,24 +130,15 @@ void SubModelView::SetFocus(const std::vector& columns_focus) { rows_focus_.push_back(i); } } - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); } -CoreModel::CoreModel(const Model* model) - : Model(*model), - full_model_(model), - core2full_row_map_(model->num_elements()), - full2core_row_map_(model->num_elements()), - core2full_col_map_(model->num_subsets()) { +CoreModel::CoreModel(const Model* model) : Model(), full_model_(model) { CHECK(ElementIndex(full_model_.num_elements()) < null_element_index) << "Max element index is reserved."; CHECK(SubsetIndex(full_model_.num_subsets()) < null_subset_index) << "Max subset index is reserved."; - - absl::c_iota(core2full_row_map_, FullElementIndex()); - absl::c_iota(full2core_row_map_, ElementIndex()); - absl::c_iota(core2full_col_map_, FullSubsetIndex()); + ResetToIdentitySubModel(); } CoreModel::CoreModel(const Model* model, @@ -161,6 +157,18 @@ CoreModel::CoreModel(const Model* model, SetFocus(columns_focus); } +void CoreModel::ResetToIdentitySubModel() { + core2full_row_map_.resize(full_model_.num_elements()); + full2core_row_map_.resize(full_model_.num_elements()); + core2full_col_map_.resize(full_model_.num_subsets()); + absl::c_iota(core2full_row_map_, FullElementIndex()); + absl::c_iota(full2core_row_map_, ElementIndex()); + absl::c_iota(core2full_col_map_, FullSubsetIndex()); + fixed_cost_ = .0; + fixed_columns_.clear(); + static_cast(*this) = Model(full_model_.base()); +} + // Note: Assumes that columns_focus covers all rows for which rows_flags is true // (i.e.: non-covered rows should be set to false to rows_flags). This property // get exploited to keep the rows in the same ordering of the original model @@ -177,7 +185,6 @@ void CoreModel::SetFocus(const std::vector& columns_focus) { // Now we can fill the new core model for (FullSubsetIndex full_j : columns_focus) { - core2full_col_map_.push_back(full_j); bool first_row = true; for (FullElementIndex full_i : full_model_.columns()[full_j]) { ElementIndex core_i = full2core_row_map_[full_i]; @@ -190,10 +197,13 @@ void CoreModel::SetFocus(const std::vector& columns_focus) { submodel.AddElementToLastSubset(core_i); } } + // Handle empty columns + if (!first_row) { + core2full_col_map_.push_back(full_j); + } } submodel.CreateSparseRowView(); - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); } @@ -273,7 +283,7 @@ Model CoreModel::MakeNewCoreModel( return new_submodel; } -Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { +Cost CoreModel::FixMoreColumns(const std::vector& columns_to_fix) { if (columns_to_fix.empty()) { return .0; } @@ -288,7 +298,6 @@ Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { // Create new model object applying the computed mappings static_cast(*this) = MakeNewCoreModel(new_c2f_row_map); - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); DCHECK(absl::c_is_sorted(core2full_col_map_)); DCHECK(absl::c_is_sorted(core2full_row_map_)); @@ -296,4 +305,14 @@ Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { return fixed_cost_ - old_fixed_cost; } +void CoreModel::ResetColumnFixing( + const std::vector& columns_to_fix, const DualState&) { + ResetToIdentitySubModel(); + std::vector core_column_to_fix; + for (FullSubsetIndex full_j : columns_to_fix) { + core_column_to_fix.push_back(static_cast(full_j)); + } + FixMoreColumns(core_column_to_fix); +} + } // namespace operations_research::scp \ No newline at end of file diff --git a/ortools/set_cover/set_cover_submodel.h b/ortools/set_cover/set_cover_submodel.h index 18dc20d29e..13b4ce93f1 100644 --- a/ortools/set_cover/set_cover_submodel.h +++ b/ortools/set_cover/set_cover_submodel.h @@ -26,6 +26,7 @@ using Model = SetCoverModel; // Forward declarations, see below for the definition of the classes. struct PrimalDualState; struct Solution; +struct DualState; // The CFT algorithm generates sub-models in two distinct ways: // @@ -112,15 +113,29 @@ class SubModelView : public IndexListModelView { // Fix the provided columns, removing them for the submodel. Rows now covered // by fixed columns are also removed from the submodel along with non-fixed // columns that only cover those rows. - virtual Cost FixColumns(const std::vector& columns_to_fix); + virtual Cost FixMoreColumns(const std::vector& columns_to_fix); + + virtual void ResetColumnFixing( + const std::vector& columns_to_fix, + const DualState& state); // Hook function for specializations. This function can be used to define a // "small" core model considering a subset of the full model through the use // of column-generation or by only selecting columns with good reduced cost in // the full model. - virtual bool UpdateCore(PrimalDualState& core_state) { return false; } + virtual bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) { + return false; + } + + StrongModelView StrongTypedFullModelView() const { + return StrongModelView(full_model_); + } private: + void ResetToIdentitySubModel(); + // Pointer to the original model const Model* full_model_; @@ -179,7 +194,6 @@ class CoreModel : private Model { ElementIndex MapFullToCoreElementIndex(FullElementIndex full_i) const { DCHECK(FullElementIndex() <= full_i && full_i < FullElementIndex(num_elements())); - DCHECK(full2core_row_map_[full_i] != null_element_index); return full2core_row_map_[full_i]; } FullSubsetIndex MapCoreToFullSubsetIndex(SubsetIndex core_j) const { @@ -210,18 +224,29 @@ class CoreModel : private Model { // Fix the provided columns, removing them for the submodel. Rows now covered // by fixed columns are also removed from the submodel along with non-fixed // columns that only cover those rows. - virtual Cost FixColumns(const std::vector& columns_to_fix); + virtual Cost FixMoreColumns(const std::vector& columns_to_fix); + + virtual void ResetColumnFixing( + const std::vector& columns_to_fix, + const DualState& state); // Hook function for specializations. This function can be used to define a // "small" core model considering a subset of the full model through the use // of column-generation or by only selecting columns with good reduced cost in // the full model. - virtual bool UpdateCore(PrimalDualState& core_state) { return false; } + virtual bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) { + return false; + } + + StrongModelView StrongTypedFullModelView() const { return full_model_; } private: void MarkNewFixingInMaps(const std::vector& columns_to_fix); CoreToFullElementMapVector MakeOrFillBothRowMaps(); Model MakeNewCoreModel(const CoreToFullElementMapVector& new_c2f_col_map); + void ResetToIdentitySubModel(); // Pointer to the original model StrongModelView full_model_; diff --git a/ortools/set_cover/set_cover_views.h b/ortools/set_cover/set_cover_views.h index d9c97d725c..12e30390eb 100644 --- a/ortools/set_cover/set_cover_views.h +++ b/ortools/set_cover/set_cover_views.h @@ -57,6 +57,10 @@ using FullElementCostVector = util_intops::StrongVector; using FullSubsetCostVector = util_intops::StrongVector; using FullElementBoolVector = util_intops::StrongVector; using FullSubsetBoolVector = util_intops::StrongVector; +using FullElementToIntVector = + util_intops::StrongVector; +using FullSubsetToIntVector = + util_intops::StrongVector; // When a sub-model is created, indicies are compacted to be consecutive and // strarting from 0 (to reduce memory usage). Core ElementIndex to original @@ -121,6 +125,7 @@ class StrongModelView { auto ElementRange() const -> util_intops::StrongIntRange { return {FullElementIndex(), FullElementIndex(num_elements())}; } + const SetCoverModel& base() const { return *model_; } private: const SetCoverModel* model_; @@ -180,6 +185,7 @@ class IndexListModelView { DCHECK(ElementIndex() <= i && i < ElementIndex(num_elements())); return (*rows_sizes_)[i]; } + const SetCoverModel& base() const { return *model_; } private: const SetCoverModel* model_; @@ -233,6 +239,10 @@ class FilterModelView { -> util_intops::FilterIndexRangeView { return {is_focus_row_}; } + bool IsFocusCol(SubsetIndex j) const { return (*is_focus_col_)[j]; } + bool IsFocusRow(ElementIndex i) const { return (*is_focus_row_)[i]; } + + const SetCoverModel& base() const { return *model_; } private: const SetCoverModel* model_; diff --git a/ortools/set_cover/views.h b/ortools/set_cover/views.h index ba4c5a46c3..49a31be0b0 100644 --- a/ortools/set_cover/views.h +++ b/ortools/set_cover/views.h @@ -19,6 +19,8 @@ #include +#include "ortools/graph/iterators.h" + // NOTE: It may be unexpected, but views provide a subscript operator that // directly accesses the underlying original container using the original // indices. This behavior is particularly relevant in the context of the Set @@ -180,6 +182,7 @@ class IndexListView { IndexListViewIterator end() const { return IndexListViewIterator(values_, indices_.end()); } + absl::Span base() const { return values_; } private: absl::Span values_; @@ -255,6 +258,7 @@ class ValueFilterView { << "Inactive value. Are you using relative indices?"; return value; } + absl::Span base() const { return values_; } private: absl::Span values_; @@ -329,6 +333,7 @@ class IndexFilterView { return IndexFilterViewIterator(values_.end(), is_active_->end(), is_active_->end()); } + absl::Span base() const { return values_; } // NOTE: uses indices of the original container, not the filtered one template @@ -374,6 +379,7 @@ class TwoLevelsView : public Lvl1ViewT { level2_type operator*() const { return level2_type(&(*iter_), active_items_); } + const Lvl1ViewT& base() const { return *this; } private: level1_iterator iter_; @@ -460,6 +466,7 @@ class TransformView { TransformViewIterator end() const { return TransformViewIterator(values_.end()); } + absl::Span base() const { return values_; } private: absl::Span values_; From edb5359fd921ba017db42c45d4797bc985fd1415 Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Tue, 13 May 2025 18:04:31 +0200 Subject: [PATCH 2/2] export from google3 --- .allstar/BUILD.bazel | 17 ++ ortools/lp_data/README.md | 0 ortools/pdlp/primal_dual_hybrid_gradient.cc | 219 ++++++++++++++---- .../pdlp/primal_dual_hybrid_gradient_test.cc | 149 +++++++++++- ortools/pdlp/solvers.proto | 13 ++ ortools/util/time_limit.h | 12 +- 6 files changed, 347 insertions(+), 63 deletions(-) create mode 100644 .allstar/BUILD.bazel create mode 100644 ortools/lp_data/README.md diff --git a/.allstar/BUILD.bazel b/.allstar/BUILD.bazel new file mode 100644 index 0000000000..c6f182d9ae --- /dev/null +++ b/.allstar/BUILD.bazel @@ -0,0 +1,17 @@ +# Copyright 2010-2025 Google LLC +# 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. + +exports_files( + glob(["**"]), + visibility = ["//ortools/open_source:__subpackages__"], +) diff --git a/ortools/lp_data/README.md b/ortools/lp_data/README.md new file mode 100644 index 0000000000..e69de29bb2 diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.cc b/ortools/pdlp/primal_dual_hybrid_gradient.cc index 2aaf59241e..d166ef7700 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient.cc +++ b/ortools/pdlp/primal_dual_hybrid_gradient.cc @@ -11,6 +11,31 @@ // See the License for the specific language governing permissions and // limitations under the License. +// We solve a QP, which we call the "original QP", by applying preprocessing +// including presolve and rescaling, which produces a new QP that we call the +// "working QP". We then solve the working QP using the Primal Dual Hybrid +// Gradient algorithm (PDHG). The optimality criteria are evaluated using the +// original QP. There are three main modules in this file: +// * The free function `PrimalDualHybridGradient()`, which is the user API, and +// is responsible for input validation that doesn't use +// ShardedQuadraticProgram, creating a `PreprocessSolver`, and calling +// `PreprocessSolver::PreprocessAndSolve()`. +// * The class `PreprocessSolver`, which is responsible for everything that +// touches the original QP, including input validation that uses +// ShardedQuadraticProgram, the preprocessing, converting solutions to the +// working QP back to solutions to the original QP, and termination checks. It +// also creates a `Solver` object and calls `Solver::Solve()`. +// * The class `Solver`, which is responsible for everything that only touches +// the working QP. It keeps a pointer to `PreprocessSolver` and calls methods +// on it when it needs access to the original QP, e.g. termination checks. +// When feasibility polishing is enabled the main solve's `Solver` object +// creates additional `Solver` objects periodically to do the feasibility +// polishing (in `Solver::TryPrimalPolishing()` and +// `Solver::TryDualPolishing()`). +// The main reason for having two separate classes `PreprocessSolver` and +// `Solver` is the fact that feasibility polishing mode uses a single +// `PreprocessSolver` object with multiple `Solver` objects. + #include "ortools/pdlp/primal_dual_hybrid_gradient.h" #include @@ -313,6 +338,7 @@ SolverResult ConstructSolverResult(VectorXd primal_solution, .solve_log = std::move(solve_log)}; } +// See comment at top of file. class PreprocessSolver { public: // Assumes that `qp` and `params` are valid. @@ -505,6 +531,7 @@ class PreprocessSolver { IterationStatsCallback iteration_stats_callback_; }; +// See comment at top of file. class Solver { public: // `preprocess_solver` should not be nullptr, and the `PreprocessSolver` @@ -556,6 +583,15 @@ class Solver { // cause infinite and NaN values. constexpr static double kDivergentMovement = 1.0e100; + // The total number of iterations in feasibility polishing is at most + // `4 * iterations_completed_ / kFeasibilityIterationFraction`. + // One factor of two is because there are both primal and dual feasibility + // polishing phases, and the other factor of two is because + // `next_feasibility_polishing_iteration` increases by a factor of 2 each + // feasibility polishing phase, so the sum of iteration limits is at most + // twice the last value. + constexpr static int kFeasibilityIterationFraction = 8; + // Attempts to solve primal and dual feasibility subproblems starting at the // average iterate, for at most `iteration_limit` iterations each. If // successful, returns a `SolverResult`, otherwise nullopt. Appends @@ -2255,6 +2291,36 @@ IterationStats WorkFromFeasibilityPolishing(const SolveLog& solve_log) { return result; } +bool TerminationReasonIsInterrupted(const TerminationReason reason) { + return reason == TERMINATION_REASON_INTERRUPTED_BY_USER; +} + +bool TerminationReasonIsWorkLimitNotInterrupted( + const TerminationReason reason) { + return reason == TERMINATION_REASON_ITERATION_LIMIT || + reason == TERMINATION_REASON_TIME_LIMIT || + reason == TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT; +} + +// Note: `TERMINATION_REASON_INTERRUPTED_BY_USER` is treated as a work limit +// (that was determined in real-time by the user). +bool TerminationReasonIsWorkLimit(const TerminationReason reason) { + return TerminationReasonIsWorkLimitNotInterrupted(reason) || + TerminationReasonIsInterrupted(reason); +} + +bool DoFeasibilityPolishingAfterLimitsReached( + const PrimalDualHybridGradientParams& params, + const TerminationReason reason) { + if (TerminationReasonIsWorkLimitNotInterrupted(reason)) { + return params.apply_feasibility_polishing_after_limits_reached(); + } + if (TerminationReasonIsInterrupted(reason)) { + return params.apply_feasibility_polishing_if_solver_is_interrupted(); + } + return false; +} + std::optional Solver::MajorIterationAndTerminationCheck( const IterationType iteration_type, const bool force_numerical_termination, const std::atomic* interrupt_solve, @@ -2272,12 +2338,12 @@ std::optional Solver::MajorIterationAndTerminationCheck( IterationStats stats = CreateSimpleIterationStats(restart); IterationStats full_work_stats = AddWorkStats(stats, work_from_feasibility_polishing); + std::optional simple_termination_reason = + CheckSimpleTerminationCriteria(params_.termination_criteria(), + full_work_stats, interrupt_solve); const bool check_termination = major_iteration_cycle % params_.termination_check_frequency() == 0 || - CheckSimpleTerminationCriteria(params_.termination_criteria(), - full_work_stats, interrupt_solve) - .has_value() || - force_numerical_termination; + simple_termination_reason.has_value() || force_numerical_termination; // We check termination on every major iteration. DCHECK(!is_major_iteration || check_termination); if (check_termination) { @@ -2304,6 +2370,19 @@ std::optional Solver::MajorIterationAndTerminationCheck( } // We've terminated. if (maybe_termination_reason.has_value()) { + if (iteration_type == IterationType::kNormal && + DoFeasibilityPolishingAfterLimitsReached( + params_, maybe_termination_reason->reason)) { + const std::optional feasibility_result = + TryFeasibilityPolishing( + iterations_completed_ / kFeasibilityIterationFraction, + interrupt_solve, solve_log); + if (feasibility_result.has_value()) { + LOG(INFO) << "Returning result from feasibility polishing after " + "limits reached"; + return *feasibility_result; + } + } IterationStats terminating_full_stats = AddWorkStats(stats, work_from_feasibility_polishing); return PickSolutionAndConstructSolverResult( @@ -2573,15 +2652,6 @@ FeasibilityPolishingDetails BuildFeasibilityPolishingDetails( return details; } -// Note: `TERMINATION_REASON_INTERRUPTED_BY_USER` is treated as a work limit -// (that was determined in real-time by the user). -bool TerminationReasonIsWorkLimit(const TerminationReason reason) { - return reason == TERMINATION_REASON_ITERATION_LIMIT || - reason == TERMINATION_REASON_TIME_LIMIT || - reason == TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT || - reason == TERMINATION_REASON_INTERRUPTED_BY_USER; -} - std::optional Solver::TryFeasibilityPolishing( const int iteration_limit, const std::atomic* interrupt_solve, SolveLog& solve_log) { @@ -2600,12 +2670,20 @@ std::optional Solver::TryFeasibilityPolishing( // polishing, it is usually increased, and an experiment (on MIPLIB2017) // shows that this test reduces the iteration count by 3-4% on average. if (!ObjectiveGapMet(optimality_criteria, first_convergence_info)) { - if (params_.verbosity_level() >= 2) { - SOLVER_LOG(&preprocess_solver_->Logger(), - "Skipping feasibility polishing because the objective gap " - "is too large."); + std::optional simple_termination_reason = + CheckSimpleTerminationCriteria(params_.termination_criteria(), + TotalWorkSoFar(solve_log), + interrupt_solve); + if (!(simple_termination_reason.has_value() && + DoFeasibilityPolishingAfterLimitsReached( + params_, simple_termination_reason->reason))) { + if (params_.verbosity_level() >= 2) { + SOLVER_LOG(&preprocess_solver_->Logger(), + "Skipping feasibility polishing because the objective gap " + "is too large."); + } + return std::nullopt; } - return std::nullopt; } if (params_.verbosity_level() >= 2) { @@ -2623,7 +2701,17 @@ std::optional Solver::TryFeasibilityPolishing( } if (TerminationReasonIsWorkLimit( primal_result.solve_log.termination_reason())) { - return std::nullopt; + // Have we also reached the overall work limit? If so, consider finishing + // the final polishing phase. + std::optional simple_termination_reason = + CheckSimpleTerminationCriteria(params_.termination_criteria(), + TotalWorkSoFar(solve_log), + interrupt_solve); + if (!(simple_termination_reason.has_value() && + DoFeasibilityPolishingAfterLimitsReached( + params_, simple_termination_reason->reason))) { + return std::nullopt; + } } else if (primal_result.solve_log.termination_reason() != TERMINATION_REASON_OPTIMAL) { // Note: `TERMINATION_REASON_PRIMAL_INFEASIBLE` could happen normally, but @@ -2651,9 +2739,29 @@ std::optional Solver::TryFeasibilityPolishing( TerminationReason_Name(dual_result.solve_log.termination_reason())); } + IterationStats full_stats = TotalWorkSoFar(solve_log); + std::optional simple_termination_reason = + CheckSimpleTerminationCriteria(params_.termination_criteria(), full_stats, + interrupt_solve); if (TerminationReasonIsWorkLimit( dual_result.solve_log.termination_reason())) { - return std::nullopt; + // Have we also reached the overall work limit? If so, consider falling out + // of the "if" test and returning the polishing solution anyway. + if (simple_termination_reason.has_value() && + DoFeasibilityPolishingAfterLimitsReached( + params_, simple_termination_reason->reason)) { + preprocess_solver_->ComputeConvergenceAndInfeasibilityFromWorkingSolution( + params_, primal_result.primal_solution, dual_result.dual_solution, + POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION, + full_stats.add_convergence_information(), nullptr); + return ConstructSolverResult( + std::move(primal_result.primal_solution), + std::move(dual_result.dual_solution), full_stats, + simple_termination_reason->reason, + POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION, solve_log); + } else { + return std::nullopt; + } } else if (dual_result.solve_log.termination_reason() != TERMINATION_REASON_OPTIMAL) { // Note: The comment in the corresponding location when checking the @@ -2665,7 +2773,6 @@ std::optional Solver::TryFeasibilityPolishing( return std::nullopt; } - IterationStats full_stats = TotalWorkSoFar(solve_log); preprocess_solver_->ComputeConvergenceAndInfeasibilityFromWorkingSolution( params_, primal_result.primal_solution, dual_result.dual_solution, POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION, @@ -2689,12 +2796,16 @@ std::optional Solver::TryFeasibilityPolishing( full_stats, preprocess_solver_->OriginalBoundNorms(), /*force_numerical_termination=*/false); - if (earned_termination.has_value()) { - return ConstructSolverResult(std::move(primal_result.primal_solution), - std::move(dual_result.dual_solution), - full_stats, earned_termination->reason, - POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION, - solve_log); + if (earned_termination.has_value() || + (simple_termination_reason.has_value() && + DoFeasibilityPolishingAfterLimitsReached( + params_, simple_termination_reason->reason))) { + return ConstructSolverResult( + std::move(primal_result.primal_solution), + std::move(dual_result.dual_solution), full_stats, + earned_termination.has_value() ? earned_termination->reason + : simple_termination_reason->reason, + POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION, solve_log); } // Note: A typical termination check would now call // `CheckSimpleTerminationCriteria`. However, there is no obvious iterate to @@ -2708,15 +2819,22 @@ std::optional Solver::TryFeasibilityPolishing( TerminationCriteria ReduceWorkLimitsByPreviousWork( TerminationCriteria criteria, const int iteration_limit, - const IterationStats& previous_work) { - criteria.set_iteration_limit(std::max( - 0, std::min(iteration_limit, criteria.iteration_limit() - - previous_work.iteration_number()))); - criteria.set_kkt_matrix_pass_limit( - std::max(0.0, criteria.kkt_matrix_pass_limit() - - previous_work.cumulative_kkt_matrix_passes())); - criteria.set_time_sec_limit(std::max( - 0.0, criteria.time_sec_limit() - previous_work.cumulative_time_sec())); + const IterationStats& previous_work, + bool apply_feasibility_polishing_after_limits_reached) { + if (apply_feasibility_polishing_after_limits_reached) { + criteria.set_iteration_limit(iteration_limit); + criteria.set_kkt_matrix_pass_limit(std::numeric_limits::infinity()); + criteria.set_time_sec_limit(std::numeric_limits::infinity()); + } else { + criteria.set_iteration_limit(std::max( + 0, std::min(iteration_limit, criteria.iteration_limit() - + previous_work.iteration_number()))); + criteria.set_kkt_matrix_pass_limit( + std::max(0.0, criteria.kkt_matrix_pass_limit() - + previous_work.cumulative_kkt_matrix_passes())); + criteria.set_time_sec_limit(std::max( + 0.0, criteria.time_sec_limit() - previous_work.cumulative_time_sec())); + } return criteria; } @@ -2725,9 +2843,13 @@ SolverResult Solver::TryPrimalPolishing( const std::atomic* interrupt_solve, SolveLog& solve_log) { PrimalDualHybridGradientParams primal_feasibility_params = params_; *primal_feasibility_params.mutable_termination_criteria() = - ReduceWorkLimitsByPreviousWork(params_.termination_criteria(), - iteration_limit, - TotalWorkSoFar(solve_log)); + ReduceWorkLimitsByPreviousWork( + params_.termination_criteria(), iteration_limit, + TotalWorkSoFar(solve_log), + params_.apply_feasibility_polishing_after_limits_reached()); + if (params_.apply_feasibility_polishing_if_solver_is_interrupted()) { + interrupt_solve = nullptr; + } // This will save the original objective after the swap. VectorXd objective; @@ -2785,9 +2907,13 @@ SolverResult Solver::TryDualPolishing(VectorXd starting_dual_solution, SolveLog& solve_log) { PrimalDualHybridGradientParams dual_feasibility_params = params_; *dual_feasibility_params.mutable_termination_criteria() = - ReduceWorkLimitsByPreviousWork(params_.termination_criteria(), - iteration_limit, - TotalWorkSoFar(solve_log)); + ReduceWorkLimitsByPreviousWork( + params_.termination_criteria(), iteration_limit, + TotalWorkSoFar(solve_log), + params_.apply_feasibility_polishing_after_limits_reached()); + if (params_.apply_feasibility_polishing_if_solver_is_interrupted()) { + interrupt_solve = nullptr; + } // These will initially contain the homogenous variable and constraint // bounds, but will contain the original variable and constraint bounds @@ -2883,14 +3009,6 @@ SolverResult Solver::Solve(const IterationType iteration_type, if (params_.use_feasibility_polishing() && iteration_type == IterationType::kNormal && iterations_completed_ >= next_feasibility_polishing_iteration) { - // The total number of iterations in feasibility polishing is at most - // `4 * iterations_completed_ / kFeasibilityIterationFraction`. - // One factor of two is because there are both primal and dual feasibility - // polishing phases, and the other factor of two is because - // `next_feasibility_polishing_iteration` increases by a factor of 2 each - // feasibility polishing phase, so the sum of iteration limits is at most - // twice the last value. - const int kFeasibilityIterationFraction = 8; const std::optional feasibility_result = TryFeasibilityPolishing( iterations_completed_ / kFeasibilityIterationFraction, @@ -2940,6 +3058,7 @@ SolverResult PrimalDualHybridGradient( std::move(iteration_stats_callback)); } +// See comment at top of file. SolverResult PrimalDualHybridGradient( QuadraticProgram qp, const PrimalDualHybridGradientParams& params, std::optional initial_solution, diff --git a/ortools/pdlp/primal_dual_hybrid_gradient_test.cc b/ortools/pdlp/primal_dual_hybrid_gradient_test.cc index 122ab5fd25..ef1c0e6b8b 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient_test.cc +++ b/ortools/pdlp/primal_dual_hybrid_gradient_test.cc @@ -13,13 +13,11 @@ #include "ortools/pdlp/primal_dual_hybrid_gradient.h" -#include #include #include #include #include #include -#include #include #include #include @@ -30,7 +28,6 @@ #include "absl/container/flat_hash_map.h" #include "absl/log/check.h" #include "absl/log/log.h" -#include "absl/status/status.h" #include "absl/status/statusor.h" #include "absl/strings/str_cat.h" #include "gtest/gtest.h" @@ -1520,7 +1517,6 @@ TEST(PrimalDualHybridGradientTest, EmptyQp) { } TEST(PrimalDualHybridGradientTest, RespectsInterrupt) { - std::atomic interrupt_solve; PrimalDualHybridGradientParams params; params.mutable_termination_criteria() ->mutable_simple_optimality_criteria() @@ -1529,7 +1525,7 @@ TEST(PrimalDualHybridGradientTest, RespectsInterrupt) { ->mutable_simple_optimality_criteria() ->set_eps_optimal_relative(0.0); - interrupt_solve.store(true); + std::atomic interrupt_solve = true; const SolverResult output = PrimalDualHybridGradient(TestLp(), params, &interrupt_solve); EXPECT_EQ(output.solve_log.termination_reason(), @@ -1537,7 +1533,6 @@ TEST(PrimalDualHybridGradientTest, RespectsInterrupt) { } TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) { - std::atomic interrupt_solve; PrimalDualHybridGradientParams params; params.mutable_termination_criteria() ->mutable_simple_optimality_criteria() @@ -1546,7 +1541,7 @@ TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) { ->mutable_simple_optimality_criteria() ->set_eps_optimal_relative(0.0); - interrupt_solve.store(false); + std::atomic interrupt_solve = false; auto callback = [&](const IterationCallbackInfo& info) { if (info.iteration_stats.iteration_number() >= 10) { interrupt_solve.store(true); @@ -1562,7 +1557,6 @@ TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) { } TEST(PrimalDualHybridGradientTest, IgnoresFalseInterrupt) { - std::atomic interrupt_solve; PrimalDualHybridGradientParams params; params.mutable_termination_criteria() ->mutable_simple_optimality_criteria() @@ -1572,7 +1566,7 @@ TEST(PrimalDualHybridGradientTest, IgnoresFalseInterrupt) { ->set_eps_optimal_relative(0.0); params.mutable_termination_criteria()->set_kkt_matrix_pass_limit(1); - interrupt_solve.store(false); + std::atomic interrupt_solve = false; const SolverResult output = PrimalDualHybridGradient(TestLp(), params, &interrupt_solve); EXPECT_EQ(output.solve_log.termination_reason(), @@ -1793,6 +1787,143 @@ TEST_F(FeasibilityPolishingPrimalTest, FeasibilityPolishingFindsValidSolution) { 1.0e-12)); } +TEST_F(FeasibilityPolishingPrimalTest, + NoPolishingAfterIterationLimitWhenPolishingAfterLimitsDisabled) { + // Feasibility polishing would solve the problem the first time it is + // attempted, which would be at iteration 100. + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_after_limits_reached(false); + params_.mutable_termination_criteria()->set_iteration_limit(50); + SolverResult output = PrimalDualHybridGradient(lp_, params_); + EXPECT_EQ(output.solve_log.termination_reason(), + TERMINATION_REASON_ITERATION_LIMIT); +} + +TEST_F(FeasibilityPolishingPrimalTest, + PolishingAfterIterationLimitWhenPolishingAfterLimitsEnabled) { + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_after_limits_reached(true); + params_.mutable_termination_criteria()->set_iteration_limit(50); + SolverResult output = PrimalDualHybridGradient(lp_, params_); + EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL); +} + +TEST_F(FeasibilityPolishingPrimalTest, + PolishingTerminatesAfterIterationLimitWhenPolishingAfterLimitsDisabled) { + // Feasibility polishing will be triggered at iteration 100. The iteration + // limit prevents primal polishing from completing. + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_after_limits_reached(false); + params_.mutable_termination_criteria()->set_iteration_limit(101); + SolverResult output = PrimalDualHybridGradient(lp_, params_); + EXPECT_EQ(output.solve_log.termination_reason(), + TERMINATION_REASON_ITERATION_LIMIT); +} + +TEST_F(FeasibilityPolishingPrimalTest, + PolishingContinuesAfterIterationLimitWhenPolishingAfterLimitsEnabled) { + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_after_limits_reached(true); + params_.mutable_termination_criteria()->set_iteration_limit(101); + SolverResult output = PrimalDualHybridGradient(lp_, params_); + EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL); +} + +TEST_F(FeasibilityPolishingPrimalTest, + PolishingStopsAfterContinuingAfterIterationLimitWhenNotOptimal) { + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_after_limits_reached(true); + params_.mutable_termination_criteria()->set_iteration_limit(101); + auto* opt_criteria = params_.mutable_termination_criteria() + ->mutable_detailed_optimality_criteria(); + opt_criteria->set_eps_optimal_primal_residual_absolute(1.0e-16); + opt_criteria->set_eps_optimal_primal_residual_relative(0.0); + opt_criteria->set_eps_optimal_dual_residual_absolute(1.0e-16); + opt_criteria->set_eps_optimal_dual_residual_relative(0.0); + opt_criteria->set_eps_optimal_objective_gap_absolute(1.0e-16); + opt_criteria->set_eps_optimal_objective_gap_relative(0.0); + SolverResult output = PrimalDualHybridGradient(lp_, params_); + EXPECT_EQ(output.solve_log.termination_reason(), + TERMINATION_REASON_ITERATION_LIMIT); + // 100 main iterations + at most 12 primal feasibility polishing iterations + // + at most 12 dual feasibility polishing iterations. + EXPECT_LE(output.solve_log.iteration_count(), 124); +} + +TEST_F(FeasibilityPolishingPrimalTest, + NoPolishingAfterInterruptWhenPolishingAfterInterruptDisabled) { + // Feasibility polishing would solve the problem the first time it is + // attempted, which would be at iteration 100. + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_if_solver_is_interrupted(false); + std::atomic interrupt_solve = false; + auto callback = [&](const IterationCallbackInfo& info) { + if (info.iteration_stats.iteration_number() >= 50) { + interrupt_solve.store(true); + } + }; + SolverResult output = + PrimalDualHybridGradient(lp_, params_, &interrupt_solve, + /*message_callback=*/nullptr, callback); + EXPECT_EQ(output.solve_log.termination_reason(), + TERMINATION_REASON_INTERRUPTED_BY_USER); +} + +TEST_F(FeasibilityPolishingPrimalTest, + PolishingAfterInterruptWhenPolishingAfterInterruptEnabled) { + // Feasibility polishing would solve the problem the first time it is + // attempted, which would be at iteration 100. + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_if_solver_is_interrupted(true); + std::atomic interrupt_solve = false; + auto callback = [&](const IterationCallbackInfo& info) { + if (info.iteration_stats.iteration_number() >= 50) { + interrupt_solve.store(true); + } + }; + SolverResult output = + PrimalDualHybridGradient(lp_, params_, &interrupt_solve, + /*message_callback=*/nullptr, callback); + EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL); +} + +TEST_F(FeasibilityPolishingPrimalTest, + PolishingTerminatesAfterInterruptWhenPolishingAfterInterruptDisabled) { + // Feasibility polishing would solve the problem the first time it is + // attempted, which would be at iteration 100. + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_if_solver_is_interrupted(false); + std::atomic interrupt_solve = false; + auto callback = [&](const IterationCallbackInfo& info) { + if (info.iteration_type == IterationType::kPrimalFeasibility) { + interrupt_solve.store(true); + } + }; + SolverResult output = + PrimalDualHybridGradient(lp_, params_, &interrupt_solve, + /*message_callback=*/nullptr, callback); + EXPECT_EQ(output.solve_log.termination_reason(), + TERMINATION_REASON_INTERRUPTED_BY_USER); +} + +TEST_F(FeasibilityPolishingPrimalTest, + PolishingContinuesAfterInterruptWhenPolishingAfterInterruptEnabled) { + // Feasibility polishing would solve the problem the first time it is + // attempted, which would be at iteration 100. + params_.set_use_feasibility_polishing(true); + params_.set_apply_feasibility_polishing_if_solver_is_interrupted(true); + std::atomic interrupt_solve = false; + auto callback = [&](const IterationCallbackInfo& info) { + if (info.iteration_type == IterationType::kPrimalFeasibility) { + interrupt_solve.store(true); + } + }; + SolverResult output = + PrimalDualHybridGradient(lp_, params_, &interrupt_solve, + /*message_callback=*/nullptr, callback); + EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL); +} + TEST_F(FeasibilityPolishingPrimalTest, FeasibilityPolishingDetailsInLog) { SolverResult output = PrimalDualHybridGradient(lp_, params_); diff --git a/ortools/pdlp/solvers.proto b/ortools/pdlp/solvers.proto index c87a5a90d1..f21869e9a0 100644 --- a/ortools/pdlp/solvers.proto +++ b/ortools/pdlp/solvers.proto @@ -480,5 +480,18 @@ message PrimalDualHybridGradientParams { // optional bool use_feasibility_polishing = 30 [default = false]; + // If true, feasibility polishing will be applied after the iteration limit, + // kkt limit, or time limit is reached. This can result in a solution that is + // closer to feasibility, at the expense of violating the limit by a moderate + // amount. + optional bool apply_feasibility_polishing_after_limits_reached = 33 + [default = false]; + + // If true, feasibility polishing will be applied after the solver is + // interrupted. This can result in a solution that is closer to feasibility, + // at the expense of not stopping as promptly when interrupted. + optional bool apply_feasibility_polishing_if_solver_is_interrupted = 34 + [default = false]; + reserved 13, 14, 15, 20, 21; } diff --git a/ortools/util/time_limit.h b/ortools/util/time_limit.h index cf8fbca432..2fd664b182 100644 --- a/ortools/util/time_limit.h +++ b/ortools/util/time_limit.h @@ -465,19 +465,23 @@ class NestedTimeLimit { class TimeLimitCheckEveryNCalls { public: TimeLimitCheckEveryNCalls(int N, TimeLimit* time_limit) - : time_limit_(time_limit), count_(0), frequency_(N) {} + : time_limit_(time_limit), frequency_(N) {} bool LimitReached() { if (count_++ == frequency_) { - if (time_limit_->LimitReached()) return true; + if (time_limit_->LimitReached()) { + stopped_ = true; + return true; + } count_ = 0; } - return false; + return stopped_; } private: TimeLimit* time_limit_; - int count_; + bool stopped_ = false; + int count_ = 0; const int frequency_; };