From 5bbdb200cbe57ac742ace4ee2ee9f200f3a6f570 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 16 Apr 2025 13:59:40 +0200 Subject: [PATCH 01/35] Laid down outer refinement tentative integration with existing code --- ortools/set_cover/set_cover_cft.cc | 161 +++++++++++++++++++++--- ortools/set_cover/set_cover_cft.h | 14 ++- ortools/set_cover/set_cover_submodel.cc | 31 ++++- ortools/set_cover/set_cover_submodel.h | 24 +++- ortools/set_cover/set_cover_views.h | 4 + 5 files changed, 212 insertions(+), 22 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 71628ed1b30..1e9d0e5516a 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -26,6 +26,7 @@ #include "ortools/base/stl_util.h" #include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_views.h" +#include "ortools/set_cover/views.h" namespace operations_research::scp { @@ -33,20 +34,22 @@ namespace operations_research::scp { ////////////////////////// COMMON DEFINITIONS ////////////////////////// //////////////////////////////////////////////////////////////////////// namespace { -class CoverCounters { + +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 +59,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 +70,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, @@ -648,7 +654,7 @@ 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'; @@ -754,6 +760,103 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { 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) { + 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::max(); + 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) { + auto [solution, dual_state] = RunThreePhase(model, best_state.solution); + + if (iter_counter == 0) { + best_state.dual_state = std::move(dual_state); + + // Arbitrary, assumes integral costs. Should be an algorithm parameter + cost_cutoff = best_state.dual_state.lower_bound() + .999; + } + + fix_fraction = std::min(1.0, fix_fraction * fix_increment); + if (solution.cost() < best_state.solution.cost()) { + best_state.solution = solution; + fix_fraction = fix_minimum; + } + if (solution.cost() <= cost_cutoff) { + break; + } + + 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); + } + } + + return best_state; +} + /////////////////////////////////////////////////////////////////////// //////////////////////// FULL TO CORE PRICING ///////////////////////// /////////////////////////////////////////////////////////////////////// @@ -870,7 +973,7 @@ 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) { @@ -893,13 +996,43 @@ 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) { +void FullToCoreModel::ResetColumnFixing( + const std::vector& columns_to_fix, + PrimalDualState& state) { + is_focus_col_.assign(num_subsets_, true); + is_focus_row_.assign(num_elements_, true); + StrongModelView typed_full_model = StrongTypedFullModelView(); + for (FullSubsetIndex full_j : columns_to_fix) { + IsFocusCol(full_j) = false; + for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { + IsFocusRow(full_i) = false; + } + } + for (FullSubsetIndex full_j : typed_full_model.SubsetRange()) { + if (!IsFocusCol(full_j)) { + continue; + } + IsFocusCol(full_j) = false; + for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { + if (IsFocusRow(full_i)) { + IsFocusCol(full_j) = true; + break; + } + } + } + ResetPricingPeriod(); + base::ResetColumnFixing(columns_to_fix, state); + DCHECK(FullToSubModelInvariantCheck()); +} + +bool FullToCoreModel::UpdateCore(PrimalDualState& core_state, + bool force_update) { + if (!force_update && --update_countdown_ > 0) { return false; } diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 5cdc99123eb..8b90964be1f 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -320,6 +320,13 @@ 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 ///////////////////////// /////////////////////////////////////////////////////////////////////// @@ -337,8 +344,11 @@ class FullToCoreModel : public SubModel { public: 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, + PrimalDualState& state) override; + bool UpdateCore(PrimalDualState& core_state, + bool force_update = false) override; void ResetPricingPeriod(); const DualState& best_dual_state() const { return best_dual_state_; } diff --git a/ortools/set_cover/set_cover_submodel.cc b/ortools/set_cover/set_cover_submodel.cc index 02a13eb2559..58a442ab6a3 100644 --- a/ortools/set_cover/set_cover_submodel.cc +++ b/ortools/set_cover/set_cover_submodel.cc @@ -14,6 +14,7 @@ #include "ortools/set_cover/set_cover_submodel.h" #include "ortools/base/stl_util.h" +#include "ortools/set_cover/set_cover_views.h" namespace operations_research::scp { @@ -59,7 +60,8 @@ SubModelView::SubModelView(const Model* model, SetFocus(columns_focus); } -Cost SubModelView::FixColumns(const std::vector& columns_to_fix) { +Cost SubModelView::FixMoreColumns( + const std::vector& columns_to_fix) { DCHECK(full_model_ != nullptr); Cost old_fixed_cost = fixed_cost_; if (columns_to_fix.empty()) { @@ -93,6 +95,19 @@ Cost SubModelView::FixColumns(const std::vector& columns_to_fix) { return fixed_cost_ - old_fixed_cost; } +void SubModelView::ResetColumnFixing( + const std::vector& columns_to_fix, + PrimalDualState& state) { + fixed_columns_ = columns_to_fix; + fixed_cost_ = .0; + for (FullSubsetIndex full_j : fixed_columns_) { + fixed_cost_ += + full_model_->subset_costs()[static_cast(full_j)]; + } + + UpdateCore(state, /*force update=*/true); +} + void SubModelView::SetFocus(const std::vector& columns_focus) { DCHECK(full_model_ != nullptr); DCHECK(!rows_sizes_.empty()); @@ -273,7 +288,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; } @@ -296,4 +311,16 @@ Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { return fixed_cost_ - old_fixed_cost; } +void CoreModel::ResetColumnFixing( + const std::vector& columns_to_fix, + PrimalDualState& state) { + fixed_columns_ = columns_to_fix; + fixed_cost_ = .0; + for (FullSubsetIndex full_j : fixed_columns_) { + fixed_cost_ += full_model_.subset_costs()[full_j]; + } + + UpdateCore(state, /*force update=*/true); +} + } // 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 18dc20d29e9..935d9dd87c3 100644 --- a/ortools/set_cover/set_cover_submodel.h +++ b/ortools/set_cover/set_cover_submodel.h @@ -112,13 +112,20 @@ 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, + PrimalDualState& 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(PrimalDualState& core_state, + bool force_update = false) { + return false; + } private: // Pointer to the original model @@ -210,13 +217,22 @@ 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, + PrimalDualState& 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(PrimalDualState& core_state, + bool force_update = false) { + return false; + } + + StrongModelView StrongTypedFullModelView() const { return full_model_; } private: void MarkNewFixingInMaps(const std::vector& columns_to_fix); diff --git a/ortools/set_cover/set_cover_views.h b/ortools/set_cover/set_cover_views.h index d9c97d725c1..3222aa8988e 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 From 92d534ce743590c536085ac4f8a9243164da35ec Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 16 Apr 2025 23:07:52 +0200 Subject: [PATCH 02/35] Fixed Refinement design, (still figuring out a clean way to integrate it) --- ortools/set_cover/set_cover_cft.cc | 131 ++++++++++++++---------- ortools/set_cover/set_cover_cft.h | 11 +- ortools/set_cover/set_cover_submodel.cc | 35 +++++-- ortools/set_cover/set_cover_submodel.h | 15 +-- 4 files changed, 116 insertions(+), 76 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 1e9d0e5516a..8dbb8432aec 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -19,12 +19,14 @@ #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" @@ -740,6 +742,10 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { heuristic_cbs.set_step_size(dual_bound_cbs.step_size()); std::cout << "\nHeuristic Phase:\n"; 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; } @@ -822,36 +828,41 @@ PrimalDualState RunCftHeuristic(SubModel& model, RunMultiplierBasedGreedy(model, best_state.dual_state); } - Cost cost_cutoff = std::numeric_limits::max(); + 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) { + 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); - - // Arbitrary, assumes integral costs. Should be an algorithm parameter - cost_cutoff = best_state.dual_state.lower_bound() + .999; } - - fix_fraction = std::min(1.0, fix_fraction * fix_increment); if (solution.cost() < best_state.solution.cost()) { best_state.solution = solution; fix_fraction = fix_minimum; } - if (solution.cost() <= cost_cutoff) { + cost_cutoff = std::max(cost_cutoff, + fixed_cost_before + dual_state.lower_bound()); + + if (best_state.solution.cost() - .999 <= 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); + model.ResetColumnFixing(cols_to_fix, best_state.dual_state); } + std::cout << "Fixed " << cols_to_fix.size() + << " new columns with cost: " << model.fixed_cost() << '\n'; + std::cout << "Mode sizes: rows " << model.num_focus_elements() << "/" + << model.num_elements() << ", columns " + << model.num_focus_subsets() << "/" << model.num_subsets() + << '\n'; } return best_state; @@ -1001,52 +1012,14 @@ Cost FullToCoreModel::FixMoreColumns( return fixed_cost; } -void FullToCoreModel::ResetColumnFixing( - const std::vector& columns_to_fix, - PrimalDualState& state) { - is_focus_col_.assign(num_subsets_, true); - is_focus_row_.assign(num_elements_, true); - StrongModelView typed_full_model = StrongTypedFullModelView(); - for (FullSubsetIndex full_j : columns_to_fix) { - IsFocusCol(full_j) = false; - for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { - IsFocusRow(full_i) = false; - } - } - for (FullSubsetIndex full_j : typed_full_model.SubsetRange()) { - if (!IsFocusCol(full_j)) { - continue; - } - IsFocusCol(full_j) = false; - for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { - if (IsFocusRow(full_i)) { - IsFocusCol(full_j) = true; - break; - } - } - } - ResetPricingPeriod(); - base::ResetColumnFixing(columns_to_fix, state); - DCHECK(FullToSubModelInvariantCheck()); -} - -bool FullToCoreModel::UpdateCore(PrimalDualState& core_state, - bool force_update) { - if (!force_update && --update_countdown_ > 0) { - return false; - } - - 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'; +std::vector FullToCoreModel::SelectNewCoreColumns( + const std::vector& forced_columns) { + FilterModelView fixing_full_model = FixingFullModelView(); - auto fixing_full_model = FixingFullModelView(); FullSubsetBoolVector selected_columns(fixing_full_model.num_subsets(), false); std::vector new_core_columns; - - // Always retain best solution in the core model - for (FullSubsetIndex full_j : core_state.solution.subsets()) { + // Always retain best solution in the core model (if possible) + for (FullSubsetIndex full_j : forced_columns) { if (IsFocusCol(full_j)) { new_core_columns.push_back(full_j); selected_columns[full_j] = true; @@ -1060,10 +1033,58 @@ bool FullToCoreModel::UpdateCore(PrimalDualState& core_state, // NOTE: unnecessary, but it keeps equivalence between SubModelView/SubModel absl::c_sort(new_core_columns); + return new_core_columns; +} + +void FullToCoreModel::ResetColumnFixing( + const std::vector& full_columns_to_fix, + const DualState& dual_state) { + is_focus_col_.assign(num_subsets_, true); + is_focus_row_.assign(num_elements_, true); + + full_dual_state_ = dual_state; + + // 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. + std::vector focus_columns = + SelectNewCoreColumns(full_columns_to_fix); + + // Create a new SubModel object from scratch and then fix columns + static_cast(*this) = SubModel(full_model_, focus_columns); + + // 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); +} + +bool FullToCoreModel::UpdateCore(PrimalDualState& core_state) { + if (--update_countdown_ > 0) { + return false; + } + + UpdateDualState(core_state.dual_state, full_dual_state_, best_dual_state_); + std::vector new_core_columns = + SelectNewCoreColumns(core_state.solution.subsets()); 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 - }); + + UpdatePricingPeriod(full_dual_state_, core_state); + std::cout << "Lower bounds: Real " << full_dual_state_.lower_bound() + << ", Core " << core_state.dual_state.lower_bound() << '\n'; + + // multipliers didn't cange, but reduced cost must be recomputed + core_state.dual_state.DualUpdate(*this, [](auto, auto) {}); DCHECK(FullToSubModelInvariantCheck()); return true; diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 8b90964be1f..146be3b2d3f 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -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(), @@ -346,9 +347,8 @@ class FullToCoreModel : public SubModel { FullToCoreModel(const Model* full_model); Cost FixMoreColumns(const std::vector& columns_to_fix) override; void ResetColumnFixing(const std::vector& columns_to_fix, - PrimalDualState& state) override; - bool UpdateCore(PrimalDualState& core_state, - bool force_update = false) override; + const DualState& state) override; + bool UpdateCore(PrimalDualState& core_state) override; void ResetPricingPeriod(); const DualState& best_dual_state() const { return best_dual_state_; } @@ -381,6 +381,9 @@ class FullToCoreModel : public SubModel { bool FullToSubModelInvariantCheck(); private: + std::vector SelectNewCoreColumns( + const std::vector& forced_columns); + const Model* full_model_; // Note: The `is_focus_col_` vector duplicates information already present in diff --git a/ortools/set_cover/set_cover_submodel.cc b/ortools/set_cover/set_cover_submodel.cc index 58a442ab6a3..526d95ffdaa 100644 --- a/ortools/set_cover/set_cover_submodel.cc +++ b/ortools/set_cover/set_cover_submodel.cc @@ -14,6 +14,7 @@ #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 { @@ -97,7 +98,7 @@ Cost SubModelView::FixMoreColumns( void SubModelView::ResetColumnFixing( const std::vector& columns_to_fix, - PrimalDualState& state) { + const DualState& state) { fixed_columns_ = columns_to_fix; fixed_cost_ = .0; for (FullSubsetIndex full_j : fixed_columns_) { @@ -105,7 +106,7 @@ void SubModelView::ResetColumnFixing( full_model_->subset_costs()[static_cast(full_j)]; } - UpdateCore(state, /*force update=*/true); + CHECK(!"Implementation not finished yet, the SetCoverModel object has not been filled"); } void SubModelView::SetFocus(const std::vector& columns_focus) { @@ -313,14 +314,34 @@ Cost CoreModel::FixMoreColumns(const std::vector& columns_to_fix) { void CoreModel::ResetColumnFixing( const std::vector& columns_to_fix, - PrimalDualState& state) { - fixed_columns_ = columns_to_fix; + const DualState& dual_state) { + core2full_col_map_.clear(); + core2full_row_map_.clear(); + full2core_row_map_.assign(full_model_.num_elements(), ElementIndex()); + fixed_cost_ = .0; - for (FullSubsetIndex full_j : fixed_columns_) { + fixed_columns_ = columns_to_fix; + for (FullSubsetIndex full_j : columns_to_fix) { fixed_cost_ += full_model_.subset_costs()[full_j]; + for (FullElementIndex full_i : full_model_.columns()[full_j]) { + full2core_row_map_[full_i] = null_element_index; + } } - - UpdateCore(state, /*force update=*/true); + for (FullElementIndex full_i : full_model_.ElementRange()) { + if (full2core_row_map_[full_i] != null_element_index) { + full2core_row_map_[full_i] = ElementIndex(core2full_row_map_.size()); + core2full_row_map_.push_back(full_i); + } + } + for (FullSubsetIndex full_j : full_model_.SubsetRange()) { + for (FullElementIndex full_i : full_model_.columns()[full_j]) { + if (full2core_row_map_[full_i] != null_element_index) { + core2full_col_map_.push_back(full_j); + break; + } + } + } + CHECK(!"Implementation not finished yet, the SetCoverModel object has not been filled"); } } // 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 935d9dd87c3..42578b0f7e7 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: // @@ -116,16 +117,13 @@ class SubModelView : public IndexListModelView { virtual void ResetColumnFixing( const std::vector& columns_to_fix, - PrimalDualState& state); + 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, - bool force_update = false) { - return false; - } + virtual bool UpdateCore(PrimalDualState& core_state) { return false; } private: // Pointer to the original model @@ -221,16 +219,13 @@ class CoreModel : private Model { virtual void ResetColumnFixing( const std::vector& columns_to_fix, - PrimalDualState& state); + 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, - bool force_update = false) { - return false; - } + virtual bool UpdateCore(PrimalDualState& core_state) { return false; } StrongModelView StrongTypedFullModelView() const { return full_model_; } From a49b8a42342f66946052d295cb65453b8eadb924 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Thu, 17 Apr 2025 12:55:36 +0200 Subject: [PATCH 03/35] Adding base object accessor function to views --- ortools/set_cover/set_cover_views.h | 4 ++++ ortools/set_cover/views.h | 7 +++++++ 2 files changed, 11 insertions(+) diff --git a/ortools/set_cover/set_cover_views.h b/ortools/set_cover/set_cover_views.h index 3222aa8988e..377c760aa7b 100644 --- a/ortools/set_cover/set_cover_views.h +++ b/ortools/set_cover/set_cover_views.h @@ -125,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_; @@ -184,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_; @@ -238,6 +240,8 @@ class FilterModelView { return {is_focus_row_}; } + const SetCoverModel& base() const { return *model_; } + private: const SetCoverModel* model_; const SubsetBoolVector* is_focus_col_; diff --git a/ortools/set_cover/views.h b/ortools/set_cover/views.h index ba4c5a46c38..49a31be0b02 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 f72daf8729b13f423aa02b353f43493914da1b4d Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Thu, 17 Apr 2025 17:40:25 +0200 Subject: [PATCH 04/35] Completed `ResetToIdentitySubModel` for `CoreModel` and `SubModelView` --- ortools/set_cover/set_cover_submodel.cc | 120 +++++++++--------------- ortools/set_cover/set_cover_submodel.h | 7 ++ 2 files changed, 51 insertions(+), 76 deletions(-) diff --git a/ortools/set_cover/set_cover_submodel.cc b/ortools/set_cover/set_cover_submodel.cc index 526d95ffdaa..cf9ef601a1e 100644 --- a/ortools/set_cover/set_cover_submodel.cc +++ b/ortools/set_cover/set_cover_submodel.cc @@ -19,34 +19,10 @@ 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)); } @@ -61,6 +37,23 @@ SubModelView::SubModelView(const Model* model, SetFocus(columns_focus); } +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); @@ -91,22 +84,18 @@ Cost SubModelView::FixMoreColumns( 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& state) { - fixed_columns_ = columns_to_fix; - fixed_cost_ = .0; - for (FullSubsetIndex full_j : fixed_columns_) { - fixed_cost_ += - full_model_->subset_costs()[static_cast(full_j)]; + 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)); } - - CHECK(!"Implementation not finished yet, the SetCoverModel object has not been filled"); + FixMoreColumns(core_column_to_fix); } void SubModelView::SetFocus(const std::vector& columns_focus) { @@ -141,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, @@ -177,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 @@ -209,7 +201,6 @@ void CoreModel::SetFocus(const std::vector& columns_focus) { } submodel.CreateSparseRowView(); - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); } @@ -304,7 +295,6 @@ Cost CoreModel::FixMoreColumns(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_)); @@ -313,35 +303,13 @@ Cost CoreModel::FixMoreColumns(const std::vector& columns_to_fix) { } void CoreModel::ResetColumnFixing( - const std::vector& columns_to_fix, - const DualState& dual_state) { - core2full_col_map_.clear(); - core2full_row_map_.clear(); - full2core_row_map_.assign(full_model_.num_elements(), ElementIndex()); - - fixed_cost_ = .0; - fixed_columns_ = columns_to_fix; + const std::vector& columns_to_fix, const DualState&) { + ResetToIdentitySubModel(); + std::vector core_column_to_fix; for (FullSubsetIndex full_j : columns_to_fix) { - fixed_cost_ += full_model_.subset_costs()[full_j]; - for (FullElementIndex full_i : full_model_.columns()[full_j]) { - full2core_row_map_[full_i] = null_element_index; - } - } - for (FullElementIndex full_i : full_model_.ElementRange()) { - if (full2core_row_map_[full_i] != null_element_index) { - full2core_row_map_[full_i] = ElementIndex(core2full_row_map_.size()); - core2full_row_map_.push_back(full_i); - } - } - for (FullSubsetIndex full_j : full_model_.SubsetRange()) { - for (FullElementIndex full_i : full_model_.columns()[full_j]) { - if (full2core_row_map_[full_i] != null_element_index) { - core2full_col_map_.push_back(full_j); - break; - } - } + core_column_to_fix.push_back(static_cast(full_j)); } - CHECK(!"Implementation not finished yet, the SetCoverModel object has not been filled"); + 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 42578b0f7e7..7a6818673ea 100644 --- a/ortools/set_cover/set_cover_submodel.h +++ b/ortools/set_cover/set_cover_submodel.h @@ -125,7 +125,13 @@ class SubModelView : public IndexListModelView { // the full model. virtual bool UpdateCore(PrimalDualState& core_state) { return false; } + StrongModelView StrongTypedFullModelView() const { + return StrongModelView(full_model_); + } + private: + void ResetToIdentitySubModel(); + // Pointer to the original model const Model* full_model_; @@ -233,6 +239,7 @@ class CoreModel : private Model { 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_; From 9a550be3ccec1ab84cdc961b94b01d0339ee9484 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Thu, 17 Apr 2025 17:42:18 +0200 Subject: [PATCH 05/35] Minor numerical issues strenghtening --- ortools/set_cover/set_cover_cft.cc | 32 +++++++++++++++++++++--------- ortools/set_cover/set_cover_cft.h | 1 + 2 files changed, 24 insertions(+), 9 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 8dbb8432aec..eee225643c5 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -23,7 +23,6 @@ #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" @@ -114,9 +113,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] + last_core_update_countdown_(0), last_min_lb_seen_(std::numeric_limits::max()), last_max_lb_seen_(.0), step_size_update_countdown_(20), // Arbitrary from [1] @@ -162,19 +162,26 @@ void BoundCBs::ComputeMultipliersDelta(const SubgradientContext& context, 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 > 1e9 && 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; + } + Cost upper_bound = context.best_solution.cost() - context.model.fixed_cost(); Cost delta = upper_bound - lower_bound; 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])); } } @@ -267,7 +274,7 @@ 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, 1e9); }); if (dual_state.lower_bound() > best_state.dual_state.lower_bound()) { best_state.dual_state = dual_state; @@ -719,10 +726,7 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { 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() << "/" @@ -737,6 +741,11 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { 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() - .999) { + break; + } + // Phase 2: search for good solutions HeuristicCBs heuristic_cbs; heuristic_cbs.set_step_size(dual_bound_cbs.step_size()); @@ -749,6 +758,10 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { 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() - .999) { + break; + } std::cout << "\n3Phase Bounds: Lower " << best_state.dual_state.lower_bound() << ", Upper " @@ -1067,6 +1080,7 @@ void FullToCoreModel::ResetColumnFixing( } DCHECK_EQ(columns_to_fix.size(), full_columns_to_fix.size()); FixMoreColumns(columns_to_fix); + DCHECK(FullToSubModelInvariantCheck()); } bool FullToCoreModel::UpdateCore(PrimalDualState& core_state) { diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 146be3b2d3f..9f39da17940 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -144,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_); From d39ebe289f0dd759399a31d1d2bdb325d6bd3049 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Sat, 19 Apr 2025 16:02:17 +0200 Subject: [PATCH 06/35] Time statistics --- ortools/set_cover/set_cover_cft.cc | 47 +++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index eee225643c5..c2ee3f4a35b 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -30,12 +30,34 @@ #include "ortools/set_cover/views.h" namespace operations_research::scp { +#define CFT_MEASURE_TIME //////////////////////////////////////////////////////////////////////// ////////////////////////// COMMON DEFINITIONS ////////////////////////// //////////////////////////////////////////////////////////////////////// namespace { +#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: @@ -248,6 +270,7 @@ 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); @@ -559,6 +582,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]; @@ -711,6 +736,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, @@ -834,6 +860,8 @@ std::vector SelectColumnByGapContribution( 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()) { @@ -876,7 +904,24 @@ PrimalDualState RunCftHeuristic(SubModel& model, << model.num_elements() << ", columns " << model.num_focus_subsets() << "/" << model.num_subsets() << '\n'; - } +#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: %.2f(%.1f%%)\n", subg_t, + 100 * subg_t / refinement_t); + printf("Greedy Heur time: %.2f(%.1f%%)\n", greedy_t, + 100 * greedy_t / refinement_t); + printf("SubG - Greedy time: %.2f(%.1f%%)\n", subg_t - greedy_t, + 100 * (subg_t - greedy_t) / refinement_t); + printf("3Phase time: %.2f(%.1f%%)\n", three_phase_t, + 100 * three_phase_t / refinement_t); + printf("3Phase - Subg time: %.2f(%.1f%%)\n", three_phase_t - subg_t, + 100 * (three_phase_t - subg_t) / refinement_t); + printf("Total CFT time: %.2f(%.1f%%)\n", refinement_t, 100.0); +#endif return best_state; } From b1ab5b3e34e22b95edd43200715cc0f97a98c336 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Mon, 21 Apr 2025 16:05:41 +0200 Subject: [PATCH 07/35] Improved subgradient stability with numerically challenging cases --- ortools/set_cover/set_cover_cft.cc | 85 ++++++++++++++++++------------ ortools/set_cover/set_cover_cft.h | 14 +++-- 2 files changed, 56 insertions(+), 43 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index c2ee3f4a35b..4aacae09961 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -30,6 +31,11 @@ #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 //////////////////////////////////////////////////////////////////////// @@ -138,7 +144,6 @@ BoundCBs::BoundCBs(const SubModel& model) exit_test_countdown_(300), // Arbitrary from [1] exit_test_period_(300), // Arbitrary from [1] step_size_(0.1), // Arbitrary from [1] - last_core_update_countdown_(0), last_min_lb_seen_(std::numeric_limits::max()), last_max_lb_seen_(.0), step_size_update_countdown_(20), // Arbitrary from [1] @@ -146,19 +151,17 @@ 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_dual_state.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_; @@ -169,28 +172,25 @@ bool BoundCBs::ExitCondition(const SubgradientContext& context) { 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_; + if (stable_direction_.empty()) { + stable_direction_ = 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]; Cost curr_i_mult = context.current_dual_state.multipliers()[i]; if ((curr_i_mult <= .0 && direction_[i] < .0) || - (curr_i_mult > 1e9 && 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]; + stable_direction_[i] = stabilization_coeff_ * direction_[i] + + (1.0 - stabilization_coeff_) * stable_direction_[i]; + + squared_norm_ += stable_direction_[i] * stable_direction_[i]; } if (squared_norm_ <= kTol) { @@ -198,16 +198,22 @@ void BoundCBs::ComputeMultipliersDelta(const SubgradientContext& context, 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_; + + auto& multipliers = context.current_dual_state.multipliers(); + auto& best_multipliers = context.best_dual_state.multipliers(); for (ElementIndex i : context.model.ElementRange()) { - delta_mults[i] = step_constant * direction_[i]; + delta_mults[i] = step_constant * stable_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); @@ -216,12 +222,8 @@ void BoundCBs::UpdateStepSize(Cost lower_bound) { Cost delta = last_max_lb_seen_ - last_min_lb_seen_; 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; - + if (gap <= .001) { // Arbitray from [1] + step_size_ *= 1.5; // Arbitray from [1] } else if (gap > .01) { // Arbitray from [1] step_size_ /= 2.0; // Arbitray from [1] } @@ -257,12 +259,12 @@ 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; + // 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; @@ -285,6 +287,14 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, .subgradient = subgradient}; size_t iter = 0; while (!cbs.ExitCondition(context)) { + // 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) { + 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()) { @@ -297,7 +307,8 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, cbs.ComputeMultipliersDelta(context, multipliers_delta); dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { - i_mult = std::clamp(i_mult + multipliers_delta[i], .0, 1e9); + 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; @@ -628,6 +639,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); } @@ -705,6 +717,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); }); @@ -768,7 +783,7 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { best_state.dual_state = curr_state.dual_state; } if (curr_state.dual_state.lower_bound() >= - best_state.solution.cost() - model.fixed_cost() - .999) { + best_state.solution.cost() - model.fixed_cost() - CFT_BOUND_EPSILON) { break; } @@ -785,7 +800,7 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { best_state.solution = curr_state.solution; } if (curr_state.dual_state.lower_bound() >= - best_state.solution.cost() - model.fixed_cost() - .999) { + best_state.solution.cost() - model.fixed_cost() - CFT_BOUND_EPSILON) { break; } @@ -887,7 +902,7 @@ PrimalDualState RunCftHeuristic(SubModel& model, cost_cutoff = std::max(cost_cutoff, fixed_cost_before + dual_state.lower_bound()); - if (best_state.solution.cost() - .999 <= cost_cutoff) { + if (best_state.solution.cost() - CFT_BOUND_EPSILON <= cost_cutoff) { break; } diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 9f39da17940..44e56f216ad 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -245,13 +245,12 @@ class BoundCBs : public SubgradientCBs { // 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 + // is set to 0.5 by default but can be adjusted. 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_; - + ElementCostVector stable_direction_; std::vector lagrangian_solution_; // Stopping condition @@ -259,10 +258,9 @@ class BoundCBs : public SubgradientCBs { BaseInt max_iter_countdown_; BaseInt exit_test_countdown_; BaseInt exit_test_period_; - BaseInt last_core_update_countdown_; // Step size - void UpdateStepSize(Cost lower_bound); + void UpdateStepSize(SubgradientContext context); Cost step_size_; Cost last_min_lb_seen_; Cost last_max_lb_seen_; From bcc34a911995a133cd2b01b47f56f342593440ef Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Mon, 21 Apr 2025 16:29:10 +0200 Subject: [PATCH 08/35] Removed subgradient stabilization (issues with core updates) --- ortools/set_cover/set_cover_cft.cc | 36 +++++++++++++++--------------- ortools/set_cover/set_cover_cft.h | 17 +------------- 2 files changed, 19 insertions(+), 34 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 4aacae09961..3636f3e0fe5 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -143,7 +143,8 @@ BoundCBs::BoundCBs(const SubModel& model) model.num_focus_elements()), // Arbitrary 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] @@ -151,7 +152,7 @@ BoundCBs::BoundCBs(const SubModel& model) {} bool BoundCBs::ExitCondition(const SubgradientContext& context) { - Cost best_lb = context.best_dual_state.lower_bound(); + 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; @@ -163,11 +164,19 @@ bool BoundCBs::ExitCondition(const SubgradientContext& context) { return true; } exit_test_countdown_ = exit_test_period_; - Cost best_lower_bound = context.best_dual_state.lower_bound(); + Cost best_lower_bound = context.best_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; + + 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 + BaseInt extension = context.model.fixed_cost() < kTol ? 4 : 1; + return unfixed_run_extension_++ >= extension; } void BoundCBs::ComputeMultipliersDelta(const SubgradientContext& context, @@ -175,22 +184,15 @@ void BoundCBs::ComputeMultipliersDelta(const SubgradientContext& context, direction_ = context.subgradient; MakeMinimalCoverageSubgradient(context, direction_); - if (stable_direction_.empty()) { - stable_direction_ = direction_; - } - squared_norm_ = .0; - // Stabilize current direction and compute squared norm for (ElementIndex i : context.model.ElementRange()) { 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; } - stable_direction_[i] = stabilization_coeff_ * direction_[i] + - (1.0 - stabilization_coeff_) * stable_direction_[i]; - squared_norm_ += stable_direction_[i] * stable_direction_[i]; + squared_norm_ += direction_[i] * direction_[i]; } if (squared_norm_ <= kTol) { @@ -202,12 +204,10 @@ void BoundCBs::ComputeMultipliersDelta(const SubgradientContext& 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_; - auto& multipliers = context.current_dual_state.multipliers(); - auto& best_multipliers = context.best_dual_state.multipliers(); for (ElementIndex i : context.model.ElementRange()) { - delta_mults[i] = step_constant * stable_direction_[i]; + delta_mults[i] = step_constant * direction_[i]; DCHECK(std::isfinite(delta_mults[i])); } } @@ -328,7 +328,7 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, << ", global " << best_state.solution.cost() << "\n"; if (cbs.UpdateCoreModel(model, best_state)) { - std::cout << "Core model has been updated.\n"; + std::cout << "Core model has been updated.\n"; dual_state = best_state.dual_state; } } @@ -1187,7 +1187,7 @@ void FullToCoreModel::UpdatePricingPeriod(const DualState& full_dual_state, void FullToCoreModel::UpdateDualState(const DualState& core_dual_state, DualState& full_dual_state, - DualState& best_dual_state) { + DualState& best_dual_state) { auto fixing_full_model = FixingFullModelView(); full_dual_state_.DualUpdate( fixing_full_model, [&](ElementIndex full_i, Cost& i_mult) { diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 44e56f216ad..286e58935d0 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -234,23 +234,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. 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 stable_direction_; std::vector lagrangian_solution_; // Stopping condition @@ -258,6 +242,7 @@ class BoundCBs : public SubgradientCBs { BaseInt max_iter_countdown_; BaseInt exit_test_countdown_; BaseInt exit_test_period_; + BaseInt unfixed_run_extension_; // Step size void UpdateStepSize(SubgradientContext context); From ec108900a80676e2de845a7d1bf7199f1475ba52 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 22 Apr 2025 16:30:11 +0200 Subject: [PATCH 09/35] Forcing core update at subgradient end --- ortools/set_cover/set_cover_cft.cc | 124 +++++++++++++++---------- ortools/set_cover/set_cover_cft.h | 33 ++++--- ortools/set_cover/set_cover_submodel.h | 12 ++- 3 files changed, 105 insertions(+), 64 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 3636f3e0fe5..8d55042070a 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -256,11 +256,11 @@ void BoundCBs::MakeMinimalCoverageSubgradient(const SubgradientContext& context, } } -bool BoundCBs::UpdateCoreModel(SubModel& core_model, - PrimalDualState& best_state) { - if (core_model.UpdateCore(best_state)) { - prev_best_lb_ = - std::min(prev_best_lb_, best_state.dual_state.lower_bound()); +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::min(prev_best_lb_, context.best_lower_bound); // 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); @@ -277,12 +277,15 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, 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; @@ -310,8 +313,9 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, 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); @@ -322,21 +326,35 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, if (++iter % 50 == 0) std::cout << "Subgradient " << iter << " -- Bounds: Lower " - << dual_state.lower_bound() << ", best " - << best_state.dual_state.lower_bound() << " - Upper " + << dual_state.lower_bound() << ", best " << best_lower_bound + << " - Upper " << best_state.solution.cost() - model.fixed_cost() << ", global " << best_state.solution.cost() << "\n"; - 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(); } } + + if (cbs.UpdateCoreModel(context, model, /*force=*/true)) { + std::cout << "Core model has been updated.\n"; + 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); + 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"; + << dual_state.lower_bound() << ", best " << best_lower_bound + << " - Upper " << best_state.solution.cost() - model.fixed_cost() + << ", global " << best_state.solution.cost() << "\n"; } //////////////////////////////////////////////////////////////////////// @@ -804,9 +822,9 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { break; } - std::cout << "\n3Phase Bounds: Lower " - << best_state.dual_state.lower_bound() << ", Upper " - << best_state.solution.cost() << '\n'; + std::cout << "\n3Phase Bounds: Residual Lower " + << curr_state.dual_state.lower_bound() + model.fixed_cost() + << ", Upper " << best_state.solution.cost() << '\n'; // Phase 3: Fix the best columns (diving) FixBestColumns(model, curr_state); @@ -814,7 +832,8 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { } std::cout << "\n\n3Phase End\n"; - std::cout << "Final Bounds: Lower " << best_state.dual_state.lower_bound() + std::cout << "Bounds: Residual Lower " + << curr_state.dual_state.lower_bound() + model.fixed_cost() << ", Upper " << best_state.solution.cost() << '\n'; return best_state; @@ -890,6 +909,7 @@ PrimalDualState RunCftHeuristic(SubModel& model, double fix_fraction = fix_minimum; for (BaseInt iter_counter = 0; model.num_elements() > 0; ++iter_counter) { + std::cout << "\n\nRefinement iteration: " << iter_counter << '\n'; Cost fixed_cost_before = model.fixed_cost(); auto [solution, dual_state] = RunThreePhase(model, best_state.solution); if (iter_counter == 0) { @@ -901,7 +921,9 @@ PrimalDualState RunCftHeuristic(SubModel& model, } cost_cutoff = std::max(cost_cutoff, fixed_cost_before + dual_state.lower_bound()); - + std::cout << "Refinement Bounds: Residual Lower " << cost_cutoff + << ", Real Lower " << best_state.dual_state.lower_bound() + << ", Upper " << best_state.solution.cost() << '\n'; if (best_state.solution.cost() - CFT_BOUND_EPSILON <= cost_cutoff) { break; } @@ -915,27 +937,29 @@ PrimalDualState RunCftHeuristic(SubModel& model, } std::cout << "Fixed " << cols_to_fix.size() << " new columns with cost: " << model.fixed_cost() << '\n'; - std::cout << "Mode sizes: rows " << model.num_focus_elements() << "/" + std::cout << "Model sizes: rows " << model.num_focus_elements() << "/" << model.num_elements() << ", columns " << model.num_focus_subsets() << "/" << model.num_subsets() << '\n'; + } + #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: %.2f(%.1f%%)\n", subg_t, + printf("Subgradient time: %8.2f (%.1f%%)\n", subg_t, 100 * subg_t / refinement_t); - printf("Greedy Heur time: %.2f(%.1f%%)\n", greedy_t, + printf("Greedy Heur time: %8.2f (%.1f%%)\n", greedy_t, 100 * greedy_t / refinement_t); - printf("SubG - Greedy time: %.2f(%.1f%%)\n", subg_t - greedy_t, + printf("SubG - Greedy time: %8.2f (%.1f%%)\n", subg_t - greedy_t, 100 * (subg_t - greedy_t) / refinement_t); - printf("3Phase time: %.2f(%.1f%%)\n", three_phase_t, + printf("3Phase time: %8.2f (%.1f%%)\n", three_phase_t, 100 * three_phase_t / refinement_t); - printf("3Phase - Subg time: %.2f(%.1f%%)\n", three_phase_t - subg_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: %.2f(%.1f%%)\n", refinement_t, 100.0); + printf("Total CFT time: %8.2f (%.1f%%)\n", refinement_t, 100.0); #endif return best_state; @@ -1143,36 +1167,36 @@ void FullToCoreModel::ResetColumnFixing( DCHECK(FullToSubModelInvariantCheck()); } -bool FullToCoreModel::UpdateCore(PrimalDualState& core_state) { - if (--update_countdown_ > 0) { +bool FullToCoreModel::UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) { + if (!force && --update_countdown_ > 0) { return false; } - UpdateDualState(core_state.dual_state, full_dual_state_, best_dual_state_); + UpdateMultipliers(best_multipliers, full_dual_state_, best_dual_state_); std::vector new_core_columns = - SelectNewCoreColumns(core_state.solution.subsets()); + SelectNewCoreColumns(best_solution.subsets()); SetFocus(new_core_columns); - UpdatePricingPeriod(full_dual_state_, core_state); - std::cout << "Lower bounds: Real " << full_dual_state_.lower_bound() - << ", Core " << core_state.dual_state.lower_bound() << '\n'; - - // multipliers didn't cange, but reduced cost must be recomputed - core_state.dual_state.DualUpdate(*this, [](auto, auto) {}); + UpdatePricingPeriod(full_dual_state_, best_lower_bound, + best_solution.cost() - fixed_cost()); + std::cout << "Core-update: Lower bounds: Real " + << full_dual_state_.lower_bound() << ", Core " << best_lower_bound + << '\n'; 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 delta = - core_state.dual_state.lower_bound() - full_dual_state.lower_bound(); - Cost ratio = DivideIfGE0(delta, core_state.solution.cost()); + 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_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) { @@ -1185,15 +1209,15 @@ void FullToCoreModel::UpdatePricingPeriod(const DualState& full_dual_state, update_countdown_ = update_period_; } -void FullToCoreModel::UpdateDualState(const DualState& core_dual_state, - DualState& full_dual_state, +void FullToCoreModel::UpdateMultipliers( + const ElementCostVector& core_multipliers, DualState& full_dual_state, DualState& best_dual_state) { 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 diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 286e58935d0..543c33bec31 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_; }; @@ -190,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; }; @@ -203,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; }; @@ -225,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, @@ -286,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; } @@ -332,7 +338,9 @@ class FullToCoreModel : public SubModel { Cost FixMoreColumns(const std::vector& columns_to_fix) override; void ResetColumnFixing(const std::vector& columns_to_fix, const DualState& state) override; - bool UpdateCore(PrimalDualState& core_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_; } @@ -344,9 +352,10 @@ 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); + void UpdateMultipliers(const ElementCostVector& core_multipliers, + DualState& full_dual_state, + DualState& best_dual_state); // Views are not composable (for now), so we can either access the full_model // with the strongly typed view or with the filtered view. @@ -366,7 +375,7 @@ class FullToCoreModel : public SubModel { private: std::vector SelectNewCoreColumns( - const std::vector& forced_columns); + const std::vector& forced_columns = {}); const Model* full_model_; diff --git a/ortools/set_cover/set_cover_submodel.h b/ortools/set_cover/set_cover_submodel.h index 7a6818673ea..e342f329345 100644 --- a/ortools/set_cover/set_cover_submodel.h +++ b/ortools/set_cover/set_cover_submodel.h @@ -123,7 +123,11 @@ class SubModelView : public IndexListModelView { // "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_); @@ -231,7 +235,11 @@ class CoreModel : private Model { // "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_; } From 13a5730e497db22e0bd435cc6211d36bd403c468 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 22 Apr 2025 17:21:18 +0200 Subject: [PATCH 10/35] Replacing std::cout with VLOGs --- ortools/set_cover/set_cover_cft.cc | 91 +++++++++++++++--------------- 1 file changed, 44 insertions(+), 47 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 8d55042070a..38dd609ea48 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -14,6 +14,7 @@ #include "ortools/set_cover/set_cover_cft.h" #include +#include #include #include #include @@ -288,8 +289,8 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, .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). @@ -324,12 +325,11 @@ 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_lower_bound - << " - Upper " - << best_state.solution.cost() - model.fixed_cost() - << ", global " << best_state.solution.cost() << "\n"; + VLOG_EVERY_POW_2(4) << "[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(context, model)) { dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { @@ -340,7 +340,6 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, } if (cbs.UpdateCoreModel(context, model, /*force=*/true)) { - std::cout << "Core model has been updated.\n"; dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { i_mult = best_multipliers[i]; }); @@ -351,10 +350,10 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, }); DCHECK_EQ(best_state.dual_state.lower_bound(), best_lower_bound); - std::cout << "Subgradient End" << " -- Bounds: Lower " - << dual_state.lower_bound() << ", best " << best_lower_bound - << " - Upper " << best_state.solution.cost() - model.fixed_cost() - << ", global " << best_state.solution.cost() << "\n"; + 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(); } //////////////////////////////////////////////////////////////////////// @@ -720,10 +719,10 @@ void FixBestColumns(SubModel& model, PrimalDualState& state) { // Fix columns and update the model 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) { @@ -778,24 +777,24 @@ 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) { ++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; @@ -808,7 +807,7 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { // 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()) { @@ -822,19 +821,19 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { break; } - std::cout << "\n3Phase Bounds: Residual Lower " - << curr_state.dual_state.lower_bound() + model.fixed_cost() - << ", 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 << "Bounds: Residual Lower " - << curr_state.dual_state.lower_bound() + model.fixed_cost() - << ", 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; } @@ -909,7 +908,7 @@ PrimalDualState RunCftHeuristic(SubModel& model, double fix_fraction = fix_minimum; for (BaseInt iter_counter = 0; model.num_elements() > 0; ++iter_counter) { - std::cout << "\n\nRefinement iteration: " << iter_counter << '\n'; + 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) { @@ -921,9 +920,9 @@ PrimalDualState RunCftHeuristic(SubModel& model, } cost_cutoff = std::max(cost_cutoff, fixed_cost_before + dual_state.lower_bound()); - std::cout << "Refinement Bounds: Residual Lower " << cost_cutoff - << ", Real Lower " << best_state.dual_state.lower_bound() - << ", Upper " << best_state.solution.cost() << '\n'; + 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; } @@ -935,12 +934,11 @@ PrimalDualState RunCftHeuristic(SubModel& model, if (!cols_to_fix.empty()) { model.ResetColumnFixing(cols_to_fix, best_state.dual_state); } - std::cout << "Fixed " << cols_to_fix.size() - << " new columns with cost: " << model.fixed_cost() << '\n'; - std::cout << "Model sizes: rows " << model.num_focus_elements() << "/" - << model.num_elements() << ", columns " - << model.num_focus_subsets() << "/" << model.num_subsets() - << '\n'; + 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 @@ -1181,9 +1179,8 @@ bool FullToCoreModel::UpdateCore(Cost best_lower_bound, UpdatePricingPeriod(full_dual_state_, best_lower_bound, best_solution.cost() - fixed_cost()); - std::cout << "Core-update: Lower bounds: Real " - << full_dual_state_.lower_bound() << ", Core " << best_lower_bound - << '\n'; + VLOG(3) << "[F2CU] Core-update: Lower bounds: Real " + << full_dual_state_.lower_bound() << ", Core " << best_lower_bound; DCHECK(FullToSubModelInvariantCheck()); return true; From a598cf83930629853f72b964ebcff01f7a9378e0 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Fri, 25 Apr 2025 11:45:29 +0200 Subject: [PATCH 11/35] Housekeeping --- ortools/set_cover/set_cover_cft.cc | 36 +++++++++++++++++------------- ortools/set_cover/set_cover_cft.h | 8 ++++--- 2 files changed, 25 insertions(+), 19 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 38dd609ea48..1db0d712968 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -15,6 +15,7 @@ #include #include +#include #include #include #include @@ -25,6 +26,7 @@ #include #include +#include "ortools/algorithms/radix_sort.h" #include "ortools/base/stl_util.h" #include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_submodel.h" @@ -165,10 +167,9 @@ bool BoundCBs::ExitCondition(const SubgradientContext& context) { return true; } exit_test_countdown_ = exit_test_period_; - Cost best_lower_bound = context.best_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; + 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; @@ -223,10 +224,12 @@ void BoundCBs::UpdateStepSize(SubgradientContext context) { Cost delta = last_max_lb_seen_ - last_min_lb_seen_; Cost gap = DivideIfGE0(delta, last_max_lb_seen_); - if (gap <= .001) { // Arbitray from [1] - step_size_ *= 1.5; // Arbitray from [1] + if (gap <= .001) { // Arbitray from [1] + step_size_ *= 1.5; // Arbitray from [1] + 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; @@ -261,7 +264,7 @@ 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::min(prev_best_lb_, context.best_lower_bound); + 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); @@ -295,6 +298,7 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, // 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; }); } @@ -325,11 +329,11 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, best_state.solution = solution; } - VLOG_EVERY_POW_2(4) << "[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(); + 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(context, model)) { dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { @@ -779,8 +783,8 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { } VLOG(2) << "[3PHS] Initial lower bound: " << best_state.dual_state.lower_bound() - << "Initial solution cost: " << best_state.solution.cost() - << "Starting 3-phase algorithm"; + << ", Initial solution cost: " << best_state.solution.cost() + << ", Starting 3-phase algorithm"; PrimalDualState curr_state = best_state; BaseInt iter_count = 0; @@ -1008,7 +1012,7 @@ void SelecteMinRedCostColumns(FilterModelView full_model, 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]; + return reduced_costs[j1] < reduced_costs[j2]; }); candidates.resize(max_size); } @@ -1024,7 +1028,7 @@ void SelecteMinRedCostColumns(FilterModelView full_model, static void SelectMinRedCostByRow(FilterModelView full_model, const SubsetCostVector& reduced_costs, std::vector& columns_map, - FullSubsetBoolVector& selected) { + FullSubsetBoolVector& selected) { DCHECK_EQ(reduced_costs.size(), full_model.num_subsets()); DCHECK_EQ(selected.size(), full_model.num_subsets()); diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 543c33bec31..223267f7db0 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -322,6 +322,9 @@ PrimalDualState RunCftHeuristic(SubModel& model, //////////////////////// 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. @@ -334,6 +337,7 @@ class FullToCoreModel : public SubModel { }; public: + FullToCoreModel() = default; FullToCoreModel(const Model* full_model); Cost FixMoreColumns(const std::vector& columns_to_fix) override; void ResetColumnFixing(const std::vector& columns_to_fix, @@ -376,6 +380,7 @@ class FullToCoreModel : public SubModel { private: std::vector SelectNewCoreColumns( const std::vector& forced_columns = {}); + void SizeUpdate(); const Model* full_model_; @@ -405,9 +410,6 @@ class FullToCoreModel : public SubModel { BaseInt update_max_period_; }; -// Coverage counter to decide the number of columns to keep in the core model. -static constexpr BaseInt kMinCov = 5; - } // namespace operations_research::scp #endif /* OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H */ From 79d159ff91271f6334b8c7a814e20e644815f6cb Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Sat, 26 Apr 2025 11:47:20 +0200 Subject: [PATCH 12/35] CoreUpdate more robust to potential full-model changes --- ortools/set_cover/set_cover_cft.cc | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 1db0d712968..bb509e3b127 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -1169,9 +1169,19 @@ void FullToCoreModel::ResetColumnFixing( DCHECK(FullToSubModelInvariantCheck()); } +void FullToCoreModel::SizeUpdate() { + num_subsets_ += full_model_->num_subsets() - num_subsets_; + is_focus_col_.resize(num_subsets_, true); +} + bool FullToCoreModel::UpdateCore(Cost best_lower_bound, const ElementCostVector& best_multipliers, const Solution& best_solution, bool force) { + SizeUpdate(); + if (num_focus_subsets() == FixingFullModelView().num_focus_subsets()) { + return false; + } + if (!force && --update_countdown_ > 0) { return false; } From 0f0b869c63a983a5bbef3c518cb925afd49284cf Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Sun, 27 Apr 2025 11:49:07 +0200 Subject: [PATCH 13/35] Skip CoreUpdate when columns are rougthly the same number of rows --- ortools/set_cover/set_cover_cft.cc | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index bb509e3b127..5745cece00a 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -973,9 +973,16 @@ PrimalDualState RunCftHeuristic(SubModel& model, 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()) { From b843c36194ab81aedfefaa25aed13dd600d88f3d Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Mon, 28 Apr 2025 11:51:48 +0200 Subject: [PATCH 14/35] Full-model does not require row-view anymore --- ortools/set_cover/set_cover_cft.cc | 101 +++++++++++++++-------------- 1 file changed, 52 insertions(+), 49 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 5745cece00a..cb08d04f515 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -1003,67 +1003,57 @@ std::vector ComputeTentativeFocus(StrongModelView full_model) { } void SelecteMinRedCostColumns(FilterModelView full_model, + const std::vector& candidates, + ElementToIntVector& row_coverage, + BaseInt& rows_left_to_cover, + BaseInt max_selection_size, 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) { + DCHECK(absl::c_is_sorted(candidates, [&](auto j1, auto j2) { return reduced_costs[j1] < reduced_costs[j2]; - }); - candidates.resize(max_size); - } + })); + + BaseInt selected_size = 0; for (SubsetIndex j : candidates) { + if (reduced_costs[j] > 0.1 || selected_size >= max_selection_size) { + break; + } FullSubsetIndex j_full = static_cast(j); if (!selected[j_full]) { selected[j_full] = true; new_core_columns.push_back(j_full); + for (ElementIndex i : full_model.columns()[j]) { + if (++row_coverage[i] == kMinCov) { + --rows_left_to_cover; + } + } } } } -static void SelectMinRedCostByRow(FilterModelView full_model, - const SubsetCostVector& reduced_costs, - std::vector& columns_map, +static void SelectMinRedCostByRow( + FilterModelView full_model, const std::vector& candidates, + ElementToIntVector& row_coverage, BaseInt& rows_left_to_cover, + BaseInt max_selection_size, 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()); + DCHECK(absl::c_is_sorted(candidates, [&](auto j1, auto j2) { + return reduced_costs[j1] < reduced_costs[j2]; + })); - 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; - } + for (SubsetIndex j : candidates) { + if (rows_left_to_cover == 0) { + break; } - - 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); + FullSubsetIndex j_full = static_cast(j); + if (!selected[j_full]) { + selected[j_full] = true; + new_core_columns.push_back(j_full); + for (ElementIndex i : full_model.columns()[j]) { + if (++row_coverage[i] == kMinCov) { + --rows_left_to_cover; + } } } } @@ -1132,10 +1122,23 @@ std::vector FullToCoreModel::SelectNewCoreColumns( } } - 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); + auto subset_range = fixing_full_model.SubsetRange(); + std::vector candidates(subset_range.begin(), subset_range.end()); + const auto& full_reduced_costs = full_dual_state_.reduced_costs(); + absl::c_sort(candidates, [&](auto j1, auto j2) { + return full_reduced_costs[j1] < full_reduced_costs[j2]; + }); + + ElementToIntVector row_coverage(fixing_full_model.num_elements(), 0); + BaseInt rows_left_to_cover = fixing_full_model.num_elements(); + BaseInt max_selection_size = kMinCov * fixing_full_model.num_elements(); + SelecteMinRedCostColumns(fixing_full_model, candidates, row_coverage, + rows_left_to_cover, max_selection_size, + full_reduced_costs, new_core_columns, + selected_columns); + SelectMinRedCostByRow(fixing_full_model, candidates, row_coverage, + rows_left_to_cover, max_selection_size, + full_reduced_costs, new_core_columns, selected_columns); // NOTE: unnecessary, but it keeps equivalence between SubModelView/SubModel absl::c_sort(new_core_columns); From 10de7b1bc57b9a7cb5efc317d66bc3a1c2a7c706 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 11:52:10 +0200 Subject: [PATCH 15/35] Added cft example --- ortools/set_cover/samples/cft.cc | 53 ++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) 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 00000000000..90a72ea6b0b --- /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 From 60d4c63483dbbb16581b0e93d6333d32e2d4cba0 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 17:53:31 +0200 Subject: [PATCH 16/35] More torough `FullToSubModelInvariantCheck` --- ortools/set_cover/set_cover_cft.cc | 24 ++++++++++++++++++++++++ ortools/set_cover/set_cover_cft.h | 6 ++---- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index cb08d04f515..9b232d12f60 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -1282,7 +1282,31 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { " not found in full model view.\n"); return false; } + + // Assumes corresponding elements have the same order in both models + auto core_it = sub_model.columns()[core_j].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 == sub_model.columns()[core_j].end()) { + break; + } + } + if (core_it != sub_model.columns()[core_j].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()) { FullElementIndex full_i = sub_model.MapCoreToFullElementIndex(core_i); if (!is_focus_row_[static_cast(full_i)]) { diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 223267f7db0..a4f62e19490 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -347,6 +347,8 @@ class FullToCoreModel : public SubModel { const Solution& best_solution, bool force) override; void ResetPricingPeriod(); const DualState& best_dual_state() const { return best_dual_state_; } + bool FullToSubModelInvariantCheck(); + void SizeUpdate(); private: decltype(auto) IsFocusCol(FullSubsetIndex j) { @@ -375,12 +377,8 @@ class FullToCoreModel : public SubModel { return StrongModelView(full_model_); } - bool FullToSubModelInvariantCheck(); - - private: std::vector SelectNewCoreColumns( const std::vector& forced_columns = {}); - void SizeUpdate(); const Model* full_model_; From 847d49e2e90ec911184a8f87630b0b492ab25b8b Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 17:54:08 +0200 Subject: [PATCH 17/35] Better management of empty columns --- ortools/set_cover/set_cover_cft.cc | 37 +++++++++++++++++-------- ortools/set_cover/set_cover_submodel.cc | 5 +++- ortools/set_cover/set_cover_submodel.h | 1 - 3 files changed, 29 insertions(+), 14 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 9b232d12f60..957caaf6dfe 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -26,7 +26,6 @@ #include #include -#include "ortools/algorithms/radix_sort.h" #include "ortools/base/stl_util.h" #include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_submodel.h" @@ -1019,10 +1018,10 @@ void SelecteMinRedCostColumns(FilterModelView full_model, if (reduced_costs[j] > 0.1 || selected_size >= max_selection_size) { break; } - FullSubsetIndex j_full = static_cast(j); - if (!selected[j_full]) { - selected[j_full] = true; - new_core_columns.push_back(j_full); + FullSubsetIndex full_j = static_cast(j); + if (!selected[full_j]) { + selected[full_j] = true; + new_core_columns.push_back(full_j); for (ElementIndex i : full_model.columns()[j]) { if (++row_coverage[i] == kMinCov) { --rows_left_to_cover; @@ -1046,10 +1045,10 @@ static void SelectMinRedCostByRow( if (rows_left_to_cover == 0) { break; } - FullSubsetIndex j_full = static_cast(j); - if (!selected[j_full]) { - selected[j_full] = true; - new_core_columns.push_back(j_full); + FullSubsetIndex full_j = static_cast(j); + if (!selected[full_j]) { + selected[full_j] = true; + new_core_columns.push_back(full_j); for (ElementIndex i : full_model.columns()[j]) { if (++row_coverage[i] == kMinCov) { --rows_left_to_cover; @@ -1082,13 +1081,18 @@ void FullToCoreModel::ResetPricingPeriod() { 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)) { @@ -1102,6 +1106,7 @@ Cost FullToCoreModel::FixMoreColumns( } } } + ResetPricingPeriod(); Cost fixed_cost = base::FixMoreColumns(columns_to_fix); DCHECK(FullToSubModelInvariantCheck()); @@ -1114,11 +1119,22 @@ std::vector FullToCoreModel::SelectNewCoreColumns( FullSubsetBoolVector selected_columns(fixing_full_model.num_subsets(), false); std::vector new_core_columns; + ElementToIntVector row_coverage(fixing_full_model.num_elements(), 0); + BaseInt rows_left_to_cover = fixing_full_model.num_elements(); + BaseInt max_selection_size = kMinCov * fixing_full_model.num_elements(); + // Always retain best solution in the core model (if possible) for (FullSubsetIndex full_j : forced_columns) { if (IsFocusCol(full_j)) { new_core_columns.push_back(full_j); selected_columns[full_j] = true; + + SubsetIndex j = static_cast(full_j); + for (ElementIndex i : fixing_full_model.columns()[j]) { + if (++row_coverage[i] == kMinCov) { + --rows_left_to_cover; + } + } } } @@ -1129,9 +1145,6 @@ std::vector FullToCoreModel::SelectNewCoreColumns( return full_reduced_costs[j1] < full_reduced_costs[j2]; }); - ElementToIntVector row_coverage(fixing_full_model.num_elements(), 0); - BaseInt rows_left_to_cover = fixing_full_model.num_elements(); - BaseInt max_selection_size = kMinCov * fixing_full_model.num_elements(); SelecteMinRedCostColumns(fixing_full_model, candidates, row_coverage, rows_left_to_cover, max_selection_size, full_reduced_costs, new_core_columns, diff --git a/ortools/set_cover/set_cover_submodel.cc b/ortools/set_cover/set_cover_submodel.cc index cf9ef601a1e..59202e9f925 100644 --- a/ortools/set_cover/set_cover_submodel.cc +++ b/ortools/set_cover/set_cover_submodel.cc @@ -185,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]; @@ -198,6 +197,10 @@ 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(); diff --git a/ortools/set_cover/set_cover_submodel.h b/ortools/set_cover/set_cover_submodel.h index e342f329345..13b4ce93f1c 100644 --- a/ortools/set_cover/set_cover_submodel.h +++ b/ortools/set_cover/set_cover_submodel.h @@ -194,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 { From bdc0c16eb8d51f51d6be1db8bd9fdaa45a362b83 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 30 Apr 2025 00:49:10 +0200 Subject: [PATCH 18/35] Improved `FullToSubModelInvariantCheck` --- ortools/set_cover/set_cover_cft.cc | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 957caaf6dfe..a81c5dca2fe 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -1296,8 +1296,20 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { 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 = sub_model.columns()[core_j].begin(); + 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; @@ -1308,11 +1320,11 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { " but it is not the same as the full model.\n"); return false; } - if (++core_it == sub_model.columns()[core_j].end()) { + if (++core_it == core_column.end()) { break; } } - if (core_it != sub_model.columns()[core_j].end()) { + 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"); @@ -1321,6 +1333,12 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { } 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, @@ -1329,6 +1347,7 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { return false; } } + for (FullElementIndex full_i : typed_full_model.ElementRange()) { if (!IsFocusRow(full_i)) { continue; From 2ba86b920d02e02ab0052604795383606bd6693c Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 30 Apr 2025 00:49:39 +0200 Subject: [PATCH 19/35] Refactored `FullToCoreModel` column selection --- ortools/set_cover/set_cover_cft.cc | 217 +++++++++++++++------------- ortools/set_cover/set_cover_cft.h | 8 +- ortools/set_cover/set_cover_views.h | 2 + 3 files changed, 121 insertions(+), 106 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index a81c5dca2fe..406dfbc5b56 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -176,14 +176,18 @@ bool BoundCBs::ExitCondition(const SubgradientContext& context) { // (Not in [1]): During the first unfixed iteration we want to converge closer // to the optimum - BaseInt extension = context.model.fixed_cost() < kTol ? 4 : 1; + 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) { direction_ = context.subgradient; - MakeMinimalCoverageSubgradient(context, direction_); + BaseInt extension = (context.model.fixed_columns().empty() ? 2 : 1); + if (unfixed_run_extension_ < extension) { + MakeMinimalCoverageSubgradient(context, direction_); + } squared_norm_ = .0; for (ElementIndex i : context.model.ElementRange()) { @@ -1001,62 +1005,100 @@ std::vector ComputeTentativeFocus(StrongModelView full_model) { return columns_focus; } -void SelecteMinRedCostColumns(FilterModelView full_model, - const std::vector& candidates, - ElementToIntVector& row_coverage, - BaseInt& rows_left_to_cover, - BaseInt max_selection_size, - const SubsetCostVector& reduced_costs, - std::vector& new_core_columns, - FullSubsetBoolVector& selected) { - DCHECK(absl::c_is_sorted(candidates, [&](auto j1, auto j2) { - return reduced_costs[j1] < reduced_costs[j2]; - })); - - BaseInt selected_size = 0; - for (SubsetIndex j : candidates) { - if (reduced_costs[j] > 0.1 || selected_size >= max_selection_size) { - break; +class ColumnSelector { + public: + ColumnSelector(FilterModelView full_model, + const std::vector& forced_columns, + const SubsetCostVector& reduced_costs) + : full_model_(full_model), + candidates_(), + row_cover_counts_(full_model.num_elements(), 0), + rows_left_to_cover_(full_model.num_elements()), + selection_(), + selected_(full_model.num_subsets(), false), + reduced_costs_(reduced_costs) { + // 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(j); + } } - FullSubsetIndex full_j = static_cast(j); - if (!selected[full_j]) { - selected[full_j] = true; - new_core_columns.push_back(full_j); - for (ElementIndex i : full_model.columns()[j]) { - if (++row_coverage[i] == kMinCov) { - --rows_left_to_cover; - } + + 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(); + }; + + bool SelectColumn(SubsetIndex j) { + if (selected_[j]) { + return false; + } + for (ElementIndex i : full_model_.columns()[j]) { + selected_[j] = true; + if (++row_cover_counts_[i] == kMinCov) { + --rows_left_to_cover_; } } + // Skip empty comlumns + if (selected_[j]) { + selection_.push_back(static_cast(j)); + } + return selected_[j]; } -} -static void SelectMinRedCostByRow( - FilterModelView full_model, const std::vector& candidates, - ElementToIntVector& row_coverage, BaseInt& rows_left_to_cover, - BaseInt max_selection_size, const SubsetCostVector& reduced_costs, - std::vector& new_core_columns, - FullSubsetBoolVector& selected) { - DCHECK(absl::c_is_sorted(candidates, [&](auto j1, auto j2) { - return reduced_costs[j1] < reduced_costs[j2]; - })); - - for (SubsetIndex j : candidates) { - if (rows_left_to_cover == 0) { - break; + void SelecteMinRedCostColumns() { + 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(*it++) ? 1 : 0; } - FullSubsetIndex full_j = static_cast(j); - if (!selected[full_j]) { - selected[full_j] = true; - new_core_columns.push_back(full_j); - for (ElementIndex i : full_model.columns()[j]) { - if (++row_coverage[i] == kMinCov) { - --rows_left_to_cover; + first_unselected_ = it; + } + + void SelectMinRedCostByRow() { + 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]) { + if (row_cover_counts_[i] >= kMinCov) { + continue; + } + selected_[j] = true; + if (++row_cover_counts_[i] == kMinCov) { + --rows_left_to_cover_; } } + if (selected_[j]) { + selection_.push_back(static_cast(j)); + } } } -} + + void SortSelection() { absl::c_sort(selection_); } + + const std::vector& selection() const { return selection_; } + + private: + FilterModelView full_model_; + std::vector candidates_; + std::vector::const_iterator first_unselected_; + ElementToIntVector row_cover_counts_; + BaseInt rows_left_to_cover_; + + std::vector selection_; + SubsetBoolVector selected_; + const SubsetCostVector& reduced_costs_; +}; + } // namespace FullToCoreModel::FullToCoreModel(const Model* full_model) @@ -1064,8 +1106,7 @@ FullToCoreModel::FullToCoreModel(const Model* 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(); @@ -1081,7 +1122,6 @@ void FullToCoreModel::ResetPricingPeriod() { 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); @@ -1115,54 +1155,19 @@ Cost FullToCoreModel::FixMoreColumns( std::vector FullToCoreModel::SelectNewCoreColumns( const std::vector& forced_columns) { - FilterModelView fixing_full_model = FixingFullModelView(); - - FullSubsetBoolVector selected_columns(fixing_full_model.num_subsets(), false); - std::vector new_core_columns; - ElementToIntVector row_coverage(fixing_full_model.num_elements(), 0); - BaseInt rows_left_to_cover = fixing_full_model.num_elements(); - BaseInt max_selection_size = kMinCov * fixing_full_model.num_elements(); - - // Always retain best solution in the core model (if possible) - for (FullSubsetIndex full_j : forced_columns) { - if (IsFocusCol(full_j)) { - new_core_columns.push_back(full_j); - selected_columns[full_j] = true; - - SubsetIndex j = static_cast(full_j); - for (ElementIndex i : fixing_full_model.columns()[j]) { - if (++row_coverage[i] == kMinCov) { - --rows_left_to_cover; - } - } - } - } - - auto subset_range = fixing_full_model.SubsetRange(); - std::vector candidates(subset_range.begin(), subset_range.end()); - const auto& full_reduced_costs = full_dual_state_.reduced_costs(); - absl::c_sort(candidates, [&](auto j1, auto j2) { - return full_reduced_costs[j1] < full_reduced_costs[j2]; - }); - - SelecteMinRedCostColumns(fixing_full_model, candidates, row_coverage, - rows_left_to_cover, max_selection_size, - full_reduced_costs, new_core_columns, - selected_columns); - SelectMinRedCostByRow(fixing_full_model, candidates, row_coverage, - rows_left_to_cover, max_selection_size, - full_reduced_costs, new_core_columns, selected_columns); - - // NOTE: unnecessary, but it keeps equivalence between SubModelView/SubModel - absl::c_sort(new_core_columns); - return new_core_columns; + ColumnSelector new_core_columns(FixingFullModelView(), forced_columns, + full_dual_state_.reduced_costs()); + new_core_columns.SelecteMinRedCostColumns(); + new_core_columns.SelectMinRedCostByRow(); + new_core_columns.SortSelection(); + return std::move(new_core_columns).selection(); } void FullToCoreModel::ResetColumnFixing( const std::vector& full_columns_to_fix, const DualState& dual_state) { - is_focus_col_.assign(num_subsets_, true); - is_focus_row_.assign(num_elements_, true); + is_focus_col_.assign(full_model_->num_subsets(), true); + is_focus_row_.assign(full_model_->num_elements(), true); full_dual_state_ = dual_state; @@ -1170,11 +1175,15 @@ void FullToCoreModel::ResetColumnFixing( // 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. - std::vector focus_columns = - SelectNewCoreColumns(full_columns_to_fix); + ColumnSelector col_selector(FixingFullModelView(), full_columns_to_fix, + full_dual_state_.reduced_costs()); + col_selector.SelecteMinRedCostColumns(); + col_selector.SelectMinRedCostByRow(); + col_selector.SortSelection(); // Create a new SubModel object from scratch and then fix columns - static_cast(*this) = SubModel(full_model_, focus_columns); + static_cast(*this) = + SubModel(full_model_, col_selector.selection()); // TODO(anyone): Improve this. It's Inefficient but hardly a botleneck and it // also avoid storing a full->core column map. @@ -1193,8 +1202,7 @@ void FullToCoreModel::ResetColumnFixing( } void FullToCoreModel::SizeUpdate() { - num_subsets_ += full_model_->num_subsets() - num_subsets_; - is_focus_col_.resize(num_subsets_, true); + is_focus_col_.resize(full_model_->num_subsets(), true); } bool FullToCoreModel::UpdateCore(Cost best_lower_bound, @@ -1204,10 +1212,13 @@ bool FullToCoreModel::UpdateCore(Cost best_lower_bound, if (num_focus_subsets() == FixingFullModelView().num_focus_subsets()) { return false; } - if (!force && --update_countdown_ > 0) { return false; } + if (best_lower_bound == prev_best_lower_bound_) { + return false; + } + prev_best_lower_bound_ = best_lower_bound; UpdateMultipliers(best_multipliers, full_dual_state_, best_dual_state_); std::vector new_core_columns = @@ -1216,8 +1227,10 @@ bool FullToCoreModel::UpdateCore(Cost best_lower_bound, 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; + 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(); DCHECK(FullToSubModelInvariantCheck()); return true; diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index a4f62e19490..691da509841 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -369,7 +369,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. @@ -397,9 +398,8 @@ 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_; diff --git a/ortools/set_cover/set_cover_views.h b/ortools/set_cover/set_cover_views.h index 377c760aa7b..12e30390ebc 100644 --- a/ortools/set_cover/set_cover_views.h +++ b/ortools/set_cover/set_cover_views.h @@ -239,6 +239,8 @@ 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_; } From a6e0f086ac740a437c20cda215468281a960f12b Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 30 Apr 2025 13:50:31 +0200 Subject: [PATCH 20/35] Minor `FullToCoreModel` cleaning --- ortools/set_cover/set_cover_cft.cc | 51 ++++++++++++++++-------------- 1 file changed, 28 insertions(+), 23 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 406dfbc5b56..31ca2e671ea 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -1005,6 +1005,22 @@ std::vector ComputeTentativeFocus(StrongModelView full_model) { return columns_focus; } +// 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: ColumnSelector(FilterModelView full_model, @@ -1038,13 +1054,12 @@ class ColumnSelector { return false; } for (ElementIndex i : full_model_.columns()[j]) { - selected_[j] = true; + selected_[j] = true; // Detect empty columns if (++row_cover_counts_[i] == kMinCov) { --rows_left_to_cover_; } } - // Skip empty comlumns - if (selected_[j]) { + if (selected_[j]) { // Skip empty comlumns selection_.push_back(static_cast(j)); } return selected_[j]; @@ -1069,13 +1084,9 @@ class ColumnSelector { continue; } for (ElementIndex i : full_model_.columns()[j]) { - if (row_cover_counts_[i] >= kMinCov) { - continue; - } - selected_[j] = true; - if (++row_cover_counts_[i] == kMinCov) { - --rows_left_to_cover_; - } + ++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)); @@ -1153,16 +1164,6 @@ Cost FullToCoreModel::FixMoreColumns( return fixed_cost; } -std::vector FullToCoreModel::SelectNewCoreColumns( - const std::vector& forced_columns) { - ColumnSelector new_core_columns(FixingFullModelView(), forced_columns, - full_dual_state_.reduced_costs()); - new_core_columns.SelecteMinRedCostColumns(); - new_core_columns.SelectMinRedCostByRow(); - new_core_columns.SortSelection(); - return std::move(new_core_columns).selection(); -} - void FullToCoreModel::ResetColumnFixing( const std::vector& full_columns_to_fix, const DualState& dual_state) { @@ -1221,9 +1222,13 @@ bool FullToCoreModel::UpdateCore(Cost best_lower_bound, prev_best_lower_bound_ = best_lower_bound; UpdateMultipliers(best_multipliers, full_dual_state_, best_dual_state_); - std::vector new_core_columns = - SelectNewCoreColumns(best_solution.subsets()); - SetFocus(new_core_columns); + ColumnSelector new_core_columns(FixingFullModelView(), + best_solution.subsets(), + full_dual_state_.reduced_costs()); + new_core_columns.SelecteMinRedCostColumns(); + new_core_columns.SelectMinRedCostByRow(); + new_core_columns.SortSelection(); + SetFocus(new_core_columns.selection()); UpdatePricingPeriod(full_dual_state_, best_lower_bound, best_solution.cost() - fixed_cost()); From f69287e53ca6076dfdbde9d9bfe036f75f546777 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Thu, 1 May 2025 18:28:43 +0200 Subject: [PATCH 21/35] Refactored `FullToCoreModel` to simplify extention with column generation --- ortools/set_cover/set_cover_cft.cc | 219 +++++++++++++---------------- ortools/set_cover/set_cover_cft.h | 52 ++++++- 2 files changed, 143 insertions(+), 128 deletions(-) diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 31ca2e671ea..88af8d42d57 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -1004,113 +1004,87 @@ std::vector ComputeTentativeFocus(StrongModelView full_model) { absl::c_sort(columns_focus); return columns_focus; } +} // namespace -// 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: - ColumnSelector(FilterModelView full_model, - const std::vector& forced_columns, - const SubsetCostVector& reduced_costs) - : full_model_(full_model), - candidates_(), - row_cover_counts_(full_model.num_elements(), 0), - rows_left_to_cover_(full_model.num_elements()), - selection_(), - selected_(full_model.num_subsets(), false), - reduced_costs_(reduced_costs) { - // 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(j); - } +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(); - }; + 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(); - bool SelectColumn(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)); + 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_; } - return selected_[j]; } + if (selected_[j]) { // Skip empty comlumns + selection_.push_back(static_cast(j)); + } + return selected_[j]; +} - void SelecteMinRedCostColumns() { - 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(*it++) ? 1 : 0; - } - first_unselected_ = it; +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 SelectMinRedCostByRow() { - 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)); - } +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)); } } - - void SortSelection() { absl::c_sort(selection_); } - - const std::vector& selection() const { return selection_; } - - private: - FilterModelView full_model_; - std::vector candidates_; - std::vector::const_iterator first_unselected_; - ElementToIntVector row_cover_counts_; - BaseInt rows_left_to_cover_; - - std::vector selection_; - SubsetBoolVector selected_; - const SubsetCostVector& reduced_costs_; -}; - -} // namespace +} FullToCoreModel::FullToCoreModel(const Model* full_model) : SubModel(full_model, ComputeTentativeFocus(StrongModelView(full_model))), @@ -1176,15 +1150,10 @@ void FullToCoreModel::ResetColumnFixing( // 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. - ColumnSelector col_selector(FixingFullModelView(), full_columns_to_fix, - full_dual_state_.reduced_costs()); - col_selector.SelecteMinRedCostColumns(); - col_selector.SelectMinRedCostByRow(); - col_selector.SortSelection(); - - // Create a new SubModel object from scratch and then fix columns - static_cast(*this) = - SubModel(full_model_, col_selector.selection()); + const auto& selection = col_selector_.ComputeNewSelection( + FixingFullModelView(), full_columns_to_fix, + full_dual_state_.reduced_costs()); + static_cast(*this) = SubModel(full_model_, selection); // TODO(anyone): Improve this. It's Inefficient but hardly a botleneck and it // also avoid storing a full->core column map. @@ -1206,13 +1175,7 @@ void FullToCoreModel::SizeUpdate() { is_focus_col_.resize(full_model_->num_subsets(), true); } -bool FullToCoreModel::UpdateCore(Cost best_lower_bound, - const ElementCostVector& best_multipliers, - const Solution& best_solution, bool force) { - SizeUpdate(); - if (num_focus_subsets() == FixingFullModelView().num_focus_subsets()) { - return false; - } +bool FullToCoreModel::IsTimeToUpdate(Cost best_lower_bound, bool force) { if (!force && --update_countdown_ > 0) { return false; } @@ -1220,23 +1183,31 @@ bool FullToCoreModel::UpdateCore(Cost best_lower_bound, return false; } prev_best_lower_bound_ = best_lower_bound; + return true; +} - UpdateMultipliers(best_multipliers, full_dual_state_, best_dual_state_); - ColumnSelector new_core_columns(FixingFullModelView(), - best_solution.subsets(), - full_dual_state_.reduced_costs()); - new_core_columns.SelecteMinRedCostColumns(); - new_core_columns.SelectMinRedCostByRow(); - new_core_columns.SortSelection(); - SetFocus(new_core_columns.selection()); - +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; } @@ -1261,9 +1232,8 @@ void FullToCoreModel::UpdatePricingPeriod(const DualState& full_dual_state, update_countdown_ = update_period_; } -void FullToCoreModel::UpdateMultipliers( - const ElementCostVector& core_multipliers, 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) { @@ -1287,6 +1257,7 @@ void FullToCoreModel::UpdateMultipliers( 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() { diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 691da509841..9481161a570 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -336,6 +336,46 @@ 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); @@ -348,9 +388,11 @@ class FullToCoreModel : public SubModel { 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)]; } @@ -359,9 +401,8 @@ class FullToCoreModel : public SubModel { } void UpdatePricingPeriod(const DualState& full_dual_state, Cost core_lower_bound, Cost core_upper_bound); - void UpdateMultipliers(const ElementCostVector& core_multipliers, - DualState& full_dual_state, - DualState& best_dual_state); + 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. @@ -381,6 +422,7 @@ class FullToCoreModel : public SubModel { std::vector SelectNewCoreColumns( const std::vector& forced_columns = {}); + private: const Model* full_model_; // Note: The `is_focus_col_` vector duplicates information already present in @@ -406,6 +448,8 @@ class FullToCoreModel : public SubModel { BaseInt update_countdown_; BaseInt update_period_; BaseInt update_max_period_; + + ColumnSelector col_selector_; // Here to avoid reallocations }; } // namespace operations_research::scp From 6e8cf50dc5d6d373759df919be7277345c470e0f Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 23 Apr 2025 11:29:14 +0200 Subject: [PATCH 22/35] BinPacking model --- ortools/algorithms/bin_packing.cc | 47 +++++++++++++++++++++++++++++++ ortools/algorithms/bin_packing.h | 46 ++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 ortools/algorithms/bin_packing.cc create mode 100644 ortools/algorithms/bin_packing.h diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc new file mode 100644 index 00000000000..b7d7b346906 --- /dev/null +++ b/ortools/algorithms/bin_packing.cc @@ -0,0 +1,47 @@ + +// Copyright 2025 Francesco Cavaliere +// 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. + +#include "ortools/algorithms/bin_packing.h" + +#include +#include +#include + +#include "ortools/util/filelineiter.h" + +namespace operations_research { + +void BinPackingModel::set_bin_capacity(Cost capacity) { + if (capacity <= 0) { + LOG(WARNING) << "Bin capacity must be positive."; + return; + } + bin_capcaity_ = capacity; +} +void BinPackingModel::AddItem(Cost weight) { + if (weight > bin_capcaity_) { + LOG(WARNING) << "Element weight exceeds bin capacity."; + return; + } + weigths_.push_back(weight); + is_sorted_ = false; +} +void BinPackingModel::SortWeights() { + if (!is_sorted_) { + absl::c_sort(weigths_); + is_sorted_ = true; + } +} + +} // namespace operations_research diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h new file mode 100644 index 00000000000..7b240f65fb5 --- /dev/null +++ b/ortools/algorithms/bin_packing.h @@ -0,0 +1,46 @@ + +// Copyright 2025 Francesco Cavaliere +// 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. + +#ifndef OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H +#define OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H + +#include +#include +#include + +#include "ortools/set_cover/base_types.h" + +namespace operations_research { + +class BinPackingModel { + public: + BinPackingModel() = default; + + Cost bin_capacity() const { return bin_capcaity_; } + void set_bin_capacity(Cost capacity); + const ElementCostVector& weights() const { return weigths_; } + void AddItem(Cost weight); + void SortWeights(); + ElementRange ItemRange() const { + return {ElementIndex(), ElementIndex(weigths_.size())}; + } + + private: + bool is_sorted_ = false; + Cost bin_capcaity_ = .0; + ElementCostVector weigths_ = {}; +}; + +} // namespace operations_research +#endif /* OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H */ From 69134d43b59af637bf5f6a810b018b6ff8572cdc Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 23 Apr 2025 11:29:35 +0200 Subject: [PATCH 23/35] Bin Packing instance readers --- ortools/algorithms/bin_packing.cc | 93 +++++++++++++++++++++++++++++++ ortools/algorithms/bin_packing.h | 3 + 2 files changed, 96 insertions(+) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index b7d7b346906..f999b2c2f11 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -44,4 +44,97 @@ void BinPackingModel::SortWeights() { } } +namespace { +template +bool SimpleAtoValue(absl::string_view str, int_type* out) { + return absl::SimpleAtoi(str, out); +} +bool SimpleAtoValue(absl::string_view str, bool* out) { + return absl::SimpleAtob(str, out); +} +bool SimpleAtoValue(absl::string_view str, double* out) { + return absl::SimpleAtod(str, out); +} +bool SimpleAtoValue(absl::string_view str, float* out) { + return absl::SimpleAtof(str, out); +} +} // namespace + +BinPackingModel ReadBpp(absl::string_view filename) { + BinPackingModel model; + BaseInt num_items = 0; + for (const std::string& line : + FileLines(filename, FileLineIterator::REMOVE_INLINE_CR | + FileLineIterator::REMOVE_BLANK_LINES)) { + if (num_items == 0) { + if (!SimpleAtoValue(line, &num_items)) { + LOG(WARNING) << "Invalid number of elements in file: " << line; + } + continue; + } + + Cost value = .0; + if (!SimpleAtoValue(line, &value)) { + LOG(WARNING) << "Invalid value in file: " << line; + continue; + } + if (model.bin_capacity() <= .0) { + model.set_bin_capacity(value); + } else { + model.AddItem(value); + } + } + DCHECK_GT(model.bin_capacity(), .0); + DCHECK_GT(model.weights().size(), .0); + DCHECK_EQ(num_items, model.weights().size()); + model.SortWeights(); + return model; +} + +BinPackingModel ReadCsp(absl::string_view filename) { + BinPackingModel model; + BaseInt num_item_types = 0; + for (const std::string& line : + FileLines(filename, FileLineIterator::REMOVE_INLINE_CR | + FileLineIterator::REMOVE_BLANK_LINES)) { + if (num_item_types == 0) { + if (!SimpleAtoValue(line, &num_item_types)) { + LOG(WARNING) << "Invalid number of elements in file: " << line; + } + continue; + } + if (model.bin_capacity() <= .0) { + Cost capacity = .0; + if (!SimpleAtoValue(line, &capacity)) { + LOG(WARNING) << "Invalid value in file: " << line; + } + model.set_bin_capacity(capacity); + continue; + } + + std::pair weight_and_demand = + absl::StrSplit(line, absl::ByAnyChar(" :\t"), absl::SkipEmpty()); + + Cost weight = .0; + if (!SimpleAtoValue(weight_and_demand.first, &weight)) { + LOG(WARNING) << "Invalid weight in file: " << line; + continue; + } + + BaseInt demand = 0; + if (!SimpleAtoValue(weight_and_demand.second, &demand)) { + LOG(WARNING) << "Invalid demand in file: " << line; + continue; + } + for (BaseInt i = 0; i < demand; ++i) { + model.AddItem(weight); + } + } + DCHECK_GT(model.bin_capacity(), .0); + DCHECK_GT(model.weights().size(), .0); + DCHECK_GE(num_item_types, model.weights().size()); + model.SortWeights(); + return model; +} + } // namespace operations_research diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 7b240f65fb5..343461d4d8e 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -42,5 +42,8 @@ class BinPackingModel { ElementCostVector weigths_ = {}; }; +BinPackingModel ReadBpp(absl::string_view filename); +BinPackingModel ReadCsp(absl::string_view filename); + } // namespace operations_research #endif /* OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H */ From a47a4270faf784065fc5a08b92697d29db982971 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 23 Apr 2025 16:27:36 +0200 Subject: [PATCH 24/35] Simple best fit greedy heuristic --- ortools/algorithms/bin_packing.cc | 24 +++++++++++++++++++++++- ortools/algorithms/bin_packing.h | 12 +++++++++++- 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index f999b2c2f11..43b98432198 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -132,9 +132,31 @@ BinPackingModel ReadCsp(absl::string_view filename) { } DCHECK_GT(model.bin_capacity(), .0); DCHECK_GT(model.weights().size(), .0); - DCHECK_GE(num_item_types, model.weights().size()); + DCHECK_GE(num_item_types, model.num_items()); model.SortWeights(); return model; } +void BestFit(const BinPackingModel& model, + const std::vector& items, PartialBins& bins_data) { + BaseInt new_bin = bins_data.bins.size(); + for (ElementIndex item : items) { + Cost item_weight = model.weights()[item]; + BaseInt selected_bin = new_bin; + for (BaseInt bin = 0; bin < bins_data.bins.size(); ++bin) { + Cost max_load = model.bin_capacity() - item_weight; + if (bins_data.loads[bin] <= max_load && + (selected_bin == new_bin || + bins_data.loads[bin] > bins_data.loads[selected_bin])) { + selected_bin = bin; + } + } + bins_data.bins[selected_bin].push_back(item); + bins_data.loads[selected_bin] += item_weight; + if (selected_bin == new_bin) { + ++new_bin; + } + } +} + } // namespace operations_research diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 343461d4d8e..59042e5767f 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -16,17 +16,20 @@ #define OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H #include +#include +#include #include #include #include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_submodel.h" namespace operations_research { class BinPackingModel { public: BinPackingModel() = default; - + BaseInt num_items() const { return weigths_.size(); } Cost bin_capacity() const { return bin_capcaity_; } void set_bin_capacity(Cost capacity); const ElementCostVector& weights() const { return weigths_; } @@ -42,8 +45,15 @@ class BinPackingModel { ElementCostVector weigths_ = {}; }; +struct PartialBins { + std::vector bins; + std::vector loads; +}; + BinPackingModel ReadBpp(absl::string_view filename); BinPackingModel ReadCsp(absl::string_view filename); +void BestFit(const BinPackingModel& model, + const std::vector& items, PartialBins& bins_data); } // namespace operations_research #endif /* OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H */ From 54cff1fcb3129a1b33a50df3d046b69f1661755f Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Wed, 23 Apr 2025 16:28:58 +0200 Subject: [PATCH 25/35] Simple initial bin generation --- ortools/algorithms/bin_packing.cc | 74 +++++++++++++++++++++++++++++++ ortools/algorithms/bin_packing.h | 34 ++++++++++++++ 2 files changed, 108 insertions(+) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index 43b98432198..58bcc596030 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -15,9 +15,18 @@ #include "ortools/algorithms/bin_packing.h" #include +#include #include #include +#include +#include + +#include "ortools/algorithms/radix_sort.h" +#include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_cft.h" +#include "ortools/set_cover/set_cover_model.h" +#include "ortools/set_cover/set_cover_submodel.h" #include "ortools/util/filelineiter.h" namespace operations_research { @@ -159,4 +168,69 @@ void BestFit(const BinPackingModel& model, } } +void InsertBinsIntoModel(PartialBins& bins_data, + BinPackingSetCoverModel& model) { + for (BaseInt i = 0; i < bins_data.bins.size(); ++i) { + if (bins_data.bins[i].size() > 0) { + absl::c_sort(bins_data.bins[i]); + model.AddBin(bins_data.bins[i]); + } + } +} + +BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, + PartialBins& best_solution, + BaseInt num_bins) { + BinPackingSetCoverModel scp_model; + PartialBins bins_data; + std::vector items(model.num_items()); + + absl::c_iota(items, ElementIndex(0)); + BestFit(model, items, bins_data); + if (best_solution.bins.empty() || + bins_data.bins.size() < best_solution.bins.size()) { + best_solution = bins_data; + } + InsertBinsIntoModel(bins_data, scp_model); + + // Largest first + if (!absl::c_is_sorted(model.weights(), std::greater<>())) { + bins_data.bins.clear(); + bins_data.loads.clear(); + absl::c_sort(items, [&](auto i1, auto i2) { + return model.weights()[i1] > model.weights()[i2]; + }); + BestFit(model, items, bins_data); + if (bins_data.bins.size() < best_solution.bins.size()) { + best_solution = bins_data; + } + InsertBinsIntoModel(bins_data, scp_model); + } + + std::mt19937 rnd(0); + while (scp_model.full_model().num_subsets() < num_bins) { + // Generate bins all containing a specific item + for (ElementIndex n : model.ItemRange()) { + absl::c_shuffle(items, rnd); + + auto n_it = absl::c_find(items, n); + std::iter_swap(n_it, items.end() - 1); + items.pop_back(); + + bins_data.bins.clear(); + bins_data.loads.clear(); + for (BaseInt j = 0; j < best_solution.bins.size(); ++j) { + bins_data.bins.push_back({n}); + bins_data.loads.push_back(model.weights()[n]); + BestFit(model, items, bins_data); + InsertBinsIntoModel(bins_data, scp_model); + } + + items.push_back(n); + } + } + + return scp_model; +} + } // namespace operations_research diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 59042e5767f..05d77b255b9 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -50,10 +50,44 @@ struct PartialBins { std::vector loads; }; +struct BinHash { + uint64_t operator()(const SparseColumn* bin) const { + return absl::HashOf(*bin); + } +}; + +using SubsetHashVector = util_intops::StrongVector; + +class BinPackingSetCoverModel { + public: + const scp::Model& full_model() const { return full_model_; } + + void AddBin(const SparseColumn& bin) { + DCHECK(absl::c_is_sorted(bin)); + auto it = bin_set_.find(&bin); + if (it == bin_set_.end()) { + full_model_.AddEmptySubset(1.0); + for (ElementIndex i : bin) { + full_model_.AddElementToLastSubset(i); + } + bin_set_.insert(&full_model_.columns().back()); + } + } + + private: + scp::Model full_model_; + absl::flat_hash_set bin_set_; +}; + BinPackingModel ReadBpp(absl::string_view filename); BinPackingModel ReadCsp(absl::string_view filename); void BestFit(const BinPackingModel& model, const std::vector& items, PartialBins& bins_data); + +BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, + PartialBins& best_solution, + BaseInt num_bins = 0); + } // namespace operations_research #endif /* OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H */ From 69721f90eda9707125120b0017639edf254a34b7 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Thu, 24 Apr 2025 14:42:19 +0200 Subject: [PATCH 26/35] Bin set --- ortools/algorithms/bin_packing.cc | 84 +++++++++++++++++++++++++------ ortools/algorithms/bin_packing.h | 62 +++++++++++++++-------- 2 files changed, 108 insertions(+), 38 deletions(-) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index 58bcc596030..790912a6a03 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -88,6 +89,7 @@ BinPackingModel ReadBpp(absl::string_view filename) { continue; } if (model.bin_capacity() <= .0) { + DCHECK_GT(value, .0); model.set_bin_capacity(value); } else { model.AddItem(value); @@ -148,26 +150,68 @@ BinPackingModel ReadCsp(absl::string_view filename) { void BestFit(const BinPackingModel& model, const std::vector& items, PartialBins& bins_data) { - BaseInt new_bin = bins_data.bins.size(); for (ElementIndex item : items) { Cost item_weight = model.weights()[item]; - BaseInt selected_bin = new_bin; + BaseInt selected_bin = bins_data.bins.size(); for (BaseInt bin = 0; bin < bins_data.bins.size(); ++bin) { Cost max_load = model.bin_capacity() - item_weight; if (bins_data.loads[bin] <= max_load && - (selected_bin == new_bin || + (selected_bin == bins_data.bins.size() || bins_data.loads[bin] > bins_data.loads[selected_bin])) { selected_bin = bin; } } + if (selected_bin == bins_data.bins.size()) { + bins_data.bins.emplace_back(); + bins_data.loads.emplace_back(); + } bins_data.bins[selected_bin].push_back(item); bins_data.loads[selected_bin] += item_weight; - if (selected_bin == new_bin) { - ++new_bin; + } +} + +const SparseColumn& BinPackingSetCoverModel::BinPackingModelGlobals::GetBin( + SubsetIndex j) const { + if (j < SubsetIndex(full_model.num_subsets())) { + return full_model.columns()[j]; + } + DCHECK(candidate_bin != nullptr); + return *candidate_bin; +} + +uint64_t BinPackingSetCoverModel::BinHash::operator()(SubsetIndex j) const { + DCHECK(globals != nullptr); + return absl::HashOf(globals->GetBin(j)); +} + +uint64_t BinPackingSetCoverModel::BinEq::operator()(SubsetIndex j1, + SubsetIndex j2) const { + DCHECK(globals != nullptr); + return globals->GetBin(j1) == globals->GetBin(j2); +} + +void BinPackingSetCoverModel::AddBin(const SparseColumn& bin) { + if (TryInsertBin(bin)) { + globals_.full_model.AddEmptySubset(1.0); + for (ElementIndex i : bin) { + globals_.full_model.AddElementToLastSubset(i); } } } +bool BinPackingSetCoverModel::TryInsertBin(const SparseColumn& bin) { + DCHECK(absl::c_is_sorted(bin)); + DCHECK(absl::c_adjacent_find(bin) == bin.end()); + DCHECK(globals_.candidate_bin == nullptr); + globals_.candidate_bin = &bin; + + SubsetIndex candidate_j(globals_.full_model.num_subsets()); + bool inserted = bin_set_.insert(candidate_j).second; + + globals_.candidate_bin = nullptr; + return inserted; +} + void InsertBinsIntoModel(PartialBins& bins_data, BinPackingSetCoverModel& model) { for (BaseInt i = 0; i < bins_data.bins.size(); ++i) { @@ -179,7 +223,6 @@ void InsertBinsIntoModel(PartialBins& bins_data, } BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, - PartialBins& best_solution, BaseInt num_bins) { BinPackingSetCoverModel scp_model; PartialBins bins_data; @@ -187,11 +230,9 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, absl::c_iota(items, ElementIndex(0)); BestFit(model, items, bins_data); - if (best_solution.bins.empty() || - bins_data.bins.size() < best_solution.bins.size()) { - best_solution = bins_data; - } InsertBinsIntoModel(bins_data, scp_model); + BaseInt solution_bin_num = bins_data.bins.size(); + VLOG(1) << "Best-fit solution: " << solution_bin_num << " bins"; // Largest first if (!absl::c_is_sorted(model.weights(), std::greater<>())) { @@ -201,9 +242,6 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, return model.weights()[i1] > model.weights()[i2]; }); BestFit(model, items, bins_data); - if (bins_data.bins.size() < best_solution.bins.size()) { - best_solution = bins_data; - } InsertBinsIntoModel(bins_data, scp_model); } @@ -211,6 +249,14 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, while (scp_model.full_model().num_subsets() < num_bins) { // Generate bins all containing a specific item for (ElementIndex n : model.ItemRange()) { + BaseInt unique_bin_num = scp_model.full_model().num_subsets(); + VLOG_EVERY_N_SEC(1, 5) + << "Generating bins: " << unique_bin_num << " / " << num_bins << " (" + << 100.0 * unique_bin_num / num_bins << "%)"; + if (scp_model.full_model().num_subsets() >= num_bins) { + break; + } + absl::c_shuffle(items, rnd); auto n_it = absl::c_find(items, n); @@ -219,17 +265,23 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, bins_data.bins.clear(); bins_data.loads.clear(); - for (BaseInt j = 0; j < best_solution.bins.size(); ++j) { + for (BaseInt j = 0; j < solution_bin_num; ++j) { bins_data.bins.push_back({n}); bins_data.loads.push_back(model.weights()[n]); - BestFit(model, items, bins_data); - InsertBinsIntoModel(bins_data, scp_model); } + BestFit(model, items, bins_data); + InsertBinsIntoModel(bins_data, scp_model); items.push_back(n); + + if (unique_bin_num == scp_model.full_model().num_subsets()) { + VLOG(1) << "No new bins generated."; + break; + } } } + scp_model.CompleteModel(); return scp_model; } diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 05d77b255b9..bfde95eac43 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -50,33 +50,52 @@ struct PartialBins { std::vector loads; }; -struct BinHash { - uint64_t operator()(const SparseColumn* bin) const { - return absl::HashOf(*bin); - } -}; - using SubsetHashVector = util_intops::StrongVector; class BinPackingSetCoverModel { + struct BinPackingModelGlobals { + // Dirty hack to avoid invalidation of pointers/references + // A pointer to this data structure is used to compute the hash of bins + // starting from their indices. + SetCoverModel full_model; + + // External bins do not have a valid index in the model, a temporary pointer + // is used insteda (even dirtier hack). + const SparseColumn* candidate_bin; + + const SparseColumn& GetBin(SubsetIndex j) const; + }; + + struct BinHash { + const BinPackingModelGlobals* globals; + uint64_t operator()(SubsetIndex j) const; + }; + + struct BinEq { + const BinPackingModelGlobals* globals; + uint64_t operator()(SubsetIndex j1, SubsetIndex j2) const; + }; + public: - const scp::Model& full_model() const { return full_model_; } - - void AddBin(const SparseColumn& bin) { - DCHECK(absl::c_is_sorted(bin)); - auto it = bin_set_.find(&bin); - if (it == bin_set_.end()) { - full_model_.AddEmptySubset(1.0); - for (ElementIndex i : bin) { - full_model_.AddElementToLastSubset(i); - } - bin_set_.insert(&full_model_.columns().back()); - } - } + BinPackingSetCoverModel() + : globals_{SetCoverModel(), nullptr}, + bin_set_({}, 0, BinHash{&globals_}, BinEq{&globals_}) {} + const SetCoverModel& full_model() const { return globals_.full_model; } + + void AddBin(const SparseColumn& bin); + void CompleteModel() { globals_.full_model.CreateSparseRowView(); } private: - scp::Model full_model_; - absl::flat_hash_set bin_set_; + bool TryInsertBin(const SparseColumn& bin); + + BinPackingModelGlobals globals_; + + // Contains bin indices, but it really should contains "bins" (aka, + // SparseColumn). However to avoid redundant allocations (in SetCoverModel and + // in the set) we cannot store them also here. We cannot also use + // iterators/pointers/references because they can be invalidated. So we store + // bin indices and do ungodly hacky shenanigans to get the bins from them. + absl::flat_hash_set bin_set_; }; BinPackingModel ReadBpp(absl::string_view filename); @@ -86,7 +105,6 @@ void BestFit(const BinPackingModel& model, const std::vector& items, PartialBins& bins_data); BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, - PartialBins& best_solution, BaseInt num_bins = 0); } // namespace operations_research From 80b0e1c072435e80273aba12882ae09d02324d25 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 12:04:35 +0200 Subject: [PATCH 27/35] Initial work on Bin Packing - Set Cover core model generator --- ortools/algorithms/bin_packing.cc | 14 +++++++++-- ortools/algorithms/bin_packing.h | 41 +++++++++++++++++++++++-------- 2 files changed, 43 insertions(+), 12 deletions(-) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index 790912a6a03..ea89f09754f 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -190,13 +190,15 @@ uint64_t BinPackingSetCoverModel::BinEq::operator()(SubsetIndex j1, return globals->GetBin(j1) == globals->GetBin(j2); } -void BinPackingSetCoverModel::AddBin(const SparseColumn& bin) { +bool BinPackingSetCoverModel::AddBin(const SparseColumn& bin) { if (TryInsertBin(bin)) { globals_.full_model.AddEmptySubset(1.0); for (ElementIndex i : bin) { globals_.full_model.AddElementToLastSubset(i); } + return true; } + return false; } bool BinPackingSetCoverModel::TryInsertBin(const SparseColumn& bin) { @@ -224,7 +226,7 @@ void InsertBinsIntoModel(PartialBins& bins_data, BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, BaseInt num_bins) { - BinPackingSetCoverModel scp_model; + BinPackingSetCoverModel scp_model(&model); PartialBins bins_data; std::vector items(model.num_items()); @@ -285,4 +287,12 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, return scp_model; } + +bool BinPackingSetCoverModel::UpdateCore( + Cost best_lower_bound, const ElementCostVector& best_multipliers, + const scp::Solution& best_solution, bool force) { + + return scp::FullToCoreModel::UpdateCore(best_lower_bound, best_multipliers, + best_solution, force); +} } // namespace operations_research diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index bfde95eac43..90e4d1e9680 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -21,8 +21,10 @@ #include #include +#include + #include "ortools/set_cover/base_types.h" -#include "ortools/set_cover/set_cover_submodel.h" +#include "ortools/set_cover/set_cover_cft.h" namespace operations_research { @@ -52,12 +54,14 @@ struct PartialBins { using SubsetHashVector = util_intops::StrongVector; -class BinPackingSetCoverModel { + + +class BinPackingSetCoverModel : public scp::FullToCoreModel { struct BinPackingModelGlobals { // Dirty hack to avoid invalidation of pointers/references // A pointer to this data structure is used to compute the hash of bins // starting from their indices. - SetCoverModel full_model; + scp::Model full_model; // External bins do not have a valid index in the model, a temporary pointer // is used insteda (even dirtier hack). @@ -77,25 +81,42 @@ class BinPackingSetCoverModel { }; public: - BinPackingSetCoverModel() - : globals_{SetCoverModel(), nullptr}, - bin_set_({}, 0, BinHash{&globals_}, BinEq{&globals_}) {} - const SetCoverModel& full_model() const { return globals_.full_model; } + BinPackingSetCoverModel(const BinPackingModel* bpp_model) + : globals_{scp::Model(), nullptr}, + bpp_model_(bpp_model), + knapsack_solver_(), + bin_set_({}, 0, BinHash{&globals_}, BinEq{&globals_}), + column_gen_countdown_(10), + column_gen_period_(10) {} + const scp::Model& full_model() const { return globals_.full_model; } + + bool AddBin(const SparseColumn& bin); + void CompleteModel() { + globals_.full_model.CreateSparseRowView(); + static_cast(*this) = + scp::FullToCoreModel(&globals_.full_model); + } - void AddBin(const SparseColumn& bin); - void CompleteModel() { globals_.full_model.CreateSparseRowView(); } + bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const scp::Solution& best_solution, bool force) override; private: bool TryInsertBin(const SparseColumn& bin); BinPackingModelGlobals globals_; + const BinPackingModel* bpp_model_; // Contains bin indices, but it really should contains "bins" (aka, - // SparseColumn). However to avoid redundant allocations (in SetCoverModel and + // SparseColumn). However to avoid redundant allocations (in scp::Model and // in the set) we cannot store them also here. We cannot also use // iterators/pointers/references because they can be invalidated. So we store // bin indices and do ungodly hacky shenanigans to get the bins from them. absl::flat_hash_set bin_set_; + + Cost prev_lower_bound_; + BaseInt column_gen_countdown_; + BaseInt column_gen_period_; }; BinPackingModel ReadBpp(absl::string_view filename); From 3f763101d41cf315b63d59ab944d921ae553e4df Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 12:06:43 +0200 Subject: [PATCH 28/35] Simple knapsack greedy heuristic --- ortools/algorithms/bin_packing.cc | 75 +++++++++++++++++++++++++++++++ ortools/algorithms/bin_packing.h | 34 ++++++++++++++ 2 files changed, 109 insertions(+) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index ea89f09754f..e50520fbc3c 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -287,6 +287,81 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, return scp_model; } +void ExpKnap::Solve(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit) { + capacity_ = capacity; + best_delta_ = .0; + exceptions_.clear(); + maximal_exceptions_.clear(); + break_solution_.assign(profits.size(), false); + items_.resize(profits.size()); + bnb_node_countdown_ = bnb_nodes_limit; + + for (ElementIndex i; i < ElementIndex(items_.size()); ++i) { + items_[i] = {std::max(1e-6, profits[i]), weights[i], i}; + } + + absl::c_sort(items_, [](Item i1, Item i2) { + return i1.profit / i1.weight > i2.profit / i2.weight; + }); + Heuristic(items_); + maximal_exceptions_.push_back(exceptions_); + exceptions_.clear(); + VLOG(5) << "[KPCG] Heuristic solution: cost " + << break_profit_sum_ + best_delta_; + +} + + +void ExpKnap::Heuristic( + const util_intops::StrongVector& items) { + exceptions_.clear(); + + break_profit_sum_ = break_weight_sum_ = .0; + break_it_ = items.begin(); + while (break_it_ < items.end() && + break_it_->weight <= capacity_ - break_weight_sum_) { + break_profit_sum_ += break_it_->profit; + break_weight_sum_ += break_it_->weight; + break_solution_[break_it_->index] = true; + ++break_it_; + } + Cost residual = capacity_ - break_weight_sum_; + + VLOG(5) << "[KPCG] Break solution: cost " << break_profit_sum_ + << ", residual " << residual; + + Cost profit_delta_ub = residual * break_it_->profit / break_it_->weight; + if (profit_delta_ub == .0) { + return; + } + + // Try filling the residual space with less efficient (maybe smaller) items + best_delta_ = .0; + for (auto it = break_it_; it < items.end(); it++) { + if (it->weight <= residual && it->profit > best_delta_) { + exceptions_ = {it->index}; + best_delta_ = it->profit; + if (best_delta_ >= profit_delta_ub) { + return; + } + } + } + + // Try removing an item and adding the break item + Cost min_weight = break_it_->weight - residual; + for (auto it = break_it_ - 1; it >= items.begin(); it--) { + Cost profit_delta = break_it_->profit - it->profit; + if (it->weight >= min_weight && profit_delta > best_delta_) { + exceptions_ = {break_it_->index, it->index}; + best_delta_ = profit_delta; + if (best_delta_ >= profit_delta_ub) { + return; + } + } + } +} bool BinPackingSetCoverModel::UpdateCore( Cost best_lower_bound, const ElementCostVector& best_multipliers, diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 90e4d1e9680..3aca6c57c68 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -55,6 +55,39 @@ struct PartialBins { using SubsetHashVector = util_intops::StrongVector; +class ExpKnap { + public: + struct Item { + Cost profit; // profit + Cost weight; // weight + ElementIndex index; + }; + using ItemIt = std::vector::const_iterator; + + void Solve(const ElementCostVector& profits, const ElementCostVector& weights, + Cost capacity, BaseInt bnb_nodes_limit); + + + void Heuristic(const util_intops::StrongVector& items); + + ElementBoolVector break_solution() const { return break_solution_; } + + std::vector> maximal_exceptions() const { + return maximal_exceptions_; + } + + private: + Cost capacity_; // capacity + util_intops::StrongVector items_; // items + ItemIt break_it_; + Cost break_profit_sum_; + Cost break_weight_sum_; + Cost best_delta_; + std::vector exceptions_; + std::vector> maximal_exceptions_; + ElementBoolVector break_solution_; + BaseInt bnb_node_countdown_; +}; class BinPackingSetCoverModel : public scp::FullToCoreModel { struct BinPackingModelGlobals { @@ -106,6 +139,7 @@ class BinPackingSetCoverModel : public scp::FullToCoreModel { BinPackingModelGlobals globals_; const BinPackingModel* bpp_model_; + ExpKnap knapsack_solver_; // Contains bin indices, but it really should contains "bins" (aka, // SparseColumn). However to avoid redundant allocations (in scp::Model and From 7f67e625f3d4f9783100bff394d271ee7fac7c4e Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 12:08:48 +0200 Subject: [PATCH 29/35] Simple Branch and bound for Knapsack (from Pisinger) --- ortools/algorithms/bin_packing.cc | 52 +++++++++++++++++++++++++++++++ ortools/algorithms/bin_packing.h | 2 ++ 2 files changed, 54 insertions(+) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index e50520fbc3c..bf11d62c9d9 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -311,8 +311,60 @@ void ExpKnap::Solve(const ElementCostVector& profits, VLOG(5) << "[KPCG] Heuristic solution: cost " << break_profit_sum_ + best_delta_; + EleBranch(.0, break_weight_sum_ - capacity, break_it_ - 1, break_it_ + 1); } +namespace { +static Cost BoundCheck(Cost best_delta, Cost profit_delta, Cost overweight, + ExpKnap::Item item) { + Cost bound = best_delta + 0.0; // + 1.0 for integral profits + return (profit_delta - bound) * item.weight - overweight * item.profit; +} +} // namespace + +// Adapted version from David Pisinger elebranch function: +// https://hjemmesider.diku.dk/~pisinger/expknap.c +bool ExpKnap::EleBranch(Cost profit_delta, Cost overweigth, ItemIt out_item, + ItemIt in_item) { + if (bnb_node_countdown_-- <= 0) { + return false; + } + bool improved = false; + + if (overweigth <= .0) { + if (profit_delta > best_delta_) { + best_delta_ = profit_delta; + improved = true; + VLOG(5) << "[KPCG] Improved best cost " + << break_profit_sum_ + best_delta_; + } + + bool maximal = true; + while (bnb_node_countdown_ > 0 && in_item < items_.end() && + BoundCheck(best_delta_, profit_delta, overweigth, *in_item) >= 0) { + exceptions_.push_back(in_item->index); + Cost next_delta = profit_delta + in_item->profit; + Cost next_oweight = overweigth + in_item->weight; + maximal &= !EleBranch(next_delta, next_oweight, out_item, ++in_item); + exceptions_.pop_back(); + } + + if (improved && maximal) { + maximal_exceptions_.push_back(exceptions_); + } + improved |= !maximal; + } else { + while (bnb_node_countdown_ > 0 && out_item >= items_.begin() && + BoundCheck(best_delta_, profit_delta, overweigth, *out_item) >= 0) { + exceptions_.push_back(out_item->index); + Cost next_delta = profit_delta - out_item->profit; + Cost next_oweight = overweigth - out_item->weight; + improved |= EleBranch(next_delta, next_oweight, --out_item, in_item); + exceptions_.pop_back(); + } + } + return improved; +} void ExpKnap::Heuristic( const util_intops::StrongVector& items) { diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 3aca6c57c68..6fd87d6db50 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -67,6 +67,8 @@ class ExpKnap { void Solve(const ElementCostVector& profits, const ElementCostVector& weights, Cost capacity, BaseInt bnb_nodes_limit); + bool EleBranch(Cost profit_sum, Cost overweight, ItemIt out_item, + ItemIt in_item); void Heuristic(const util_intops::StrongVector& items); From 49472840559a006930527aba56fe882be2f395d5 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 12:09:10 +0200 Subject: [PATCH 30/35] Dyamic column generation in core model update --- ortools/algorithms/bin_packing.cc | 32 +++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index bf11d62c9d9..440bbeaeb95 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -418,6 +418,38 @@ void ExpKnap::Heuristic( bool BinPackingSetCoverModel::UpdateCore( Cost best_lower_bound, const ElementCostVector& best_multipliers, const scp::Solution& best_solution, bool force) { + if (--column_gen_countdown_ <= 0 && best_lower_bound != prev_lower_bound_) { + column_gen_countdown_ = column_gen_period_; + prev_lower_bound_ = best_lower_bound; + knapsack_solver_.Solve(best_multipliers, bpp_model_->weights(), + bpp_model_->bin_capacity(), + /*bnb_nodes_limit=*/1024); + const auto& exception_list = knapsack_solver_.maximal_exceptions(); + ElementBoolVector solution = knapsack_solver_.break_solution(); + BaseInt num_added_bins = 0; + SparseColumn bin; + for (const std::vector& exception : exception_list) { + for (ElementIndex i : exception) { + solution[i] = !solution[i]; + } + + bin.clear(); + for (ElementIndex i : globals_.full_model.ElementRange()) { + if (solution[i]) { + bin.push_back(i); + } + } + num_added_bins += AddBin(bin) ? 1 : 0; + + for (ElementIndex i : exception) { + solution[i] = !solution[i]; + } + } + if (num_added_bins > 0) { + VLOG(4) << "[KPCG] Added " << num_added_bins << " / " + << globals_.full_model.num_subsets() << " bins"; + } + } return scp::FullToCoreModel::UpdateCore(best_lower_bound, best_multipliers, best_solution, force); From 9dea4093085d9ca438f5c23a05a6982dc76a79eb Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 12:09:35 +0200 Subject: [PATCH 31/35] CFT based bin packing heuristic example --- ortools/algorithms/samples/bin_packing_cft.cc | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 ortools/algorithms/samples/bin_packing_cft.cc diff --git a/ortools/algorithms/samples/bin_packing_cft.cc b/ortools/algorithms/samples/bin_packing_cft.cc new file mode 100644 index 00000000000..51195fea791 --- /dev/null +++ b/ortools/algorithms/samples/bin_packing_cft.cc @@ -0,0 +1,47 @@ +#include + +#include + +#include "ortools/algorithms/bin_packing.h" +#include "ortools/base/init_google.h" +#include "ortools/set_cover/set_cover_cft.h" + +using namespace operations_research; +ABSL_FLAG(std::string, instance, "", "BPP instance int RAIL format."); +ABSL_FLAG(int, bins, 1000, "Number of bins to generate."); + +int main(int argc, char** argv) { + InitGoogle(argv[0], &argc, &argv, true); + + BinPackingModel model = ReadBpp(absl::GetFlag(FLAGS_instance)); + + BinPackingSetCoverModel scp_model = + GenerateBins(model, absl::GetFlag(FLAGS_bins)); + + scp::PrimalDualState result = scp::RunCftHeuristic(scp_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'; + } + + if (dual.multipliers().empty()) { + std::cerr << "Error: failed to find any dual\n"; + } else { + std::cout << "Core Lower bound: " << dual.lower_bound() << '\n'; + } + + // The lower bound computed on the full model is not a real lower bound unless + // the knapsack subproblem failed to fined any negative reduced cost bin to + // add to the set cover model. + // TODO(anyone): add a flag to indicate if a valid LB has been found or not. + if (scp_model.best_dual_state().multipliers().empty()) { + std::cerr << "Error: no real dual state has been computed\n"; + } else { + std::cout << "Restricted Lower bound: " + << scp_model.best_dual_state().lower_bound() << '\n'; + } + + return EXIT_SUCCESS; +} \ No newline at end of file From a225bba05d84c3204ebb362114b5e8a1cb0a416c Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 17:46:27 +0200 Subject: [PATCH 32/35] Split initial bin generation with randomized one --- ortools/algorithms/bin_packing.cc | 55 +++++++++++++++++-------------- ortools/algorithms/bin_packing.h | 7 ++-- 2 files changed, 35 insertions(+), 27 deletions(-) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index 440bbeaeb95..b898c95811d 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -23,10 +23,8 @@ #include #include -#include "ortools/algorithms/radix_sort.h" #include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_cft.h" -#include "ortools/set_cover/set_cover_model.h" #include "ortools/set_cover/set_cover_submodel.h" #include "ortools/util/filelineiter.h" @@ -224,30 +222,12 @@ void InsertBinsIntoModel(PartialBins& bins_data, } } -BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, - BaseInt num_bins) { - BinPackingSetCoverModel scp_model(&model); +void AddRandomizedBins(const BinPackingModel& model, BaseInt num_bins, + BinPackingSetCoverModel& scp_model, std::mt19937& rnd) { PartialBins bins_data; std::vector items(model.num_items()); - absl::c_iota(items, ElementIndex(0)); - BestFit(model, items, bins_data); - InsertBinsIntoModel(bins_data, scp_model); - BaseInt solution_bin_num = bins_data.bins.size(); - VLOG(1) << "Best-fit solution: " << solution_bin_num << " bins"; - // Largest first - if (!absl::c_is_sorted(model.weights(), std::greater<>())) { - bins_data.bins.clear(); - bins_data.loads.clear(); - absl::c_sort(items, [&](auto i1, auto i2) { - return model.weights()[i1] > model.weights()[i2]; - }); - BestFit(model, items, bins_data); - InsertBinsIntoModel(bins_data, scp_model); - } - - std::mt19937 rnd(0); while (scp_model.full_model().num_subsets() < num_bins) { // Generate bins all containing a specific item for (ElementIndex n : model.ItemRange()) { @@ -267,7 +247,7 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, bins_data.bins.clear(); bins_data.loads.clear(); - for (BaseInt j = 0; j < solution_bin_num; ++j) { + for (BaseInt j = 0; j < 10; ++j) { bins_data.bins.push_back({n}); bins_data.loads.push_back(model.weights()[n]); } @@ -283,6 +263,31 @@ BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, } } + scp_model.CompleteModel(); +} + +BinPackingSetCoverModel GenerateInitialBins(const BinPackingModel& model) { + BinPackingSetCoverModel scp_model(&model); + PartialBins bins_data; + std::vector items(model.num_items()); + + absl::c_iota(items, ElementIndex(0)); + BestFit(model, items, bins_data); + InsertBinsIntoModel(bins_data, scp_model); + BaseInt solution_bin_num = bins_data.bins.size(); + VLOG(1) << "Best-fit solution: " << solution_bin_num << " bins"; + + // Largest first + if (!absl::c_is_sorted(model.weights(), std::greater<>())) { + bins_data.bins.clear(); + bins_data.loads.clear(); + absl::c_sort(items, [&](auto i1, auto i2) { + return model.weights()[i1] > model.weights()[i2]; + }); + BestFit(model, items, bins_data); + InsertBinsIntoModel(bins_data, scp_model); + } + scp_model.CompleteModel(); return scp_model; } @@ -423,7 +428,7 @@ bool BinPackingSetCoverModel::UpdateCore( prev_lower_bound_ = best_lower_bound; knapsack_solver_.Solve(best_multipliers, bpp_model_->weights(), bpp_model_->bin_capacity(), - /*bnb_nodes_limit=*/1024); + /*bnb_nodes_limit=*/10000); const auto& exception_list = knapsack_solver_.maximal_exceptions(); ElementBoolVector solution = knapsack_solver_.break_solution(); BaseInt num_added_bins = 0; @@ -451,6 +456,8 @@ bool BinPackingSetCoverModel::UpdateCore( } } + scp::FullToCoreModel::SizeUpdate(); + scp::FullToCoreModel::FullToSubModelInvariantCheck(); return scp::FullToCoreModel::UpdateCore(best_lower_bound, best_multipliers, best_solution, force); } diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 6fd87d6db50..0e627b6f5e0 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -54,7 +54,6 @@ struct PartialBins { using SubsetHashVector = util_intops::StrongVector; - class ExpKnap { public: struct Item { @@ -161,8 +160,10 @@ BinPackingModel ReadCsp(absl::string_view filename); void BestFit(const BinPackingModel& model, const std::vector& items, PartialBins& bins_data); -BinPackingSetCoverModel GenerateBins(const BinPackingModel& model, - BaseInt num_bins = 0); +BinPackingSetCoverModel GenerateInitialBins(const BinPackingModel& model); + +void AddRandomizedBins(const BinPackingModel& model, BaseInt num_bins, + BinPackingSetCoverModel& scp_model, std::mt19937& rnd); } // namespace operations_research #endif /* OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H */ From d74962d55c0ed5c98103fca29761d889169c0521 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 29 Apr 2025 17:47:05 +0200 Subject: [PATCH 33/35] Updated `bin_packing_cft.cc` example --- ortools/algorithms/samples/bin_packing_cft.cc | 32 ++++++++++++++++--- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/ortools/algorithms/samples/bin_packing_cft.cc b/ortools/algorithms/samples/bin_packing_cft.cc index 51195fea791..2b297c6cf6b 100644 --- a/ortools/algorithms/samples/bin_packing_cft.cc +++ b/ortools/algorithms/samples/bin_packing_cft.cc @@ -1,6 +1,21 @@ + +// Copyright 2025 Francesco Cavaliere +// 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. + #include #include +#include #include "ortools/algorithms/bin_packing.h" #include "ortools/base/init_google.h" @@ -15,11 +30,20 @@ int main(int argc, char** argv) { BinPackingModel model = ReadBpp(absl::GetFlag(FLAGS_instance)); - BinPackingSetCoverModel scp_model = - GenerateBins(model, absl::GetFlag(FLAGS_bins)); + // Quick run with a minimal set of bins + BinPackingSetCoverModel scp_model = GenerateInitialBins(model); + scp::PrimalDualState best_result = scp::RunCftHeuristic(scp_model); + + // Run the CFT again with more bins to get a better solution + std::mt19937 rnd(0); + AddRandomizedBins(model, absl::GetFlag(FLAGS_bins), scp_model, rnd); + scp::PrimalDualState result = + scp::RunCftHeuristic(scp_model, best_result.solution); + if (result.solution.cost() < best_result.solution.cost()) { + best_result = result; + } - scp::PrimalDualState result = scp::RunCftHeuristic(scp_model); - auto [solution, dual] = result; + auto [solution, dual] = best_result; if (solution.subsets().empty()) { std::cerr << "Error: failed to find any solution\n"; } else { From 5529e2a6b2dc0a7999273a4e7db778e264ffe30f Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Thu, 1 May 2025 18:31:19 +0200 Subject: [PATCH 34/35] Fixed some problems with BPP column generation --- ortools/algorithms/bin_packing.cc | 191 ++++++++++++++++++------------ ortools/algorithms/bin_packing.h | 31 +++-- 2 files changed, 135 insertions(+), 87 deletions(-) diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc index b898c95811d..0f512b131b9 100644 --- a/ortools/algorithms/bin_packing.cc +++ b/ortools/algorithms/bin_packing.cc @@ -21,8 +21,10 @@ #include #include +#include #include +#include "ortools/base/stl_util.h" #include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_cft.h" #include "ortools/set_cover/set_cover_submodel.h" @@ -146,13 +148,13 @@ BinPackingModel ReadCsp(absl::string_view filename) { return model; } -void BestFit(const BinPackingModel& model, +void BestFit(const ElementCostVector& weights, Cost bin_capacity, const std::vector& items, PartialBins& bins_data) { for (ElementIndex item : items) { - Cost item_weight = model.weights()[item]; + Cost item_weight = weights[item]; BaseInt selected_bin = bins_data.bins.size(); for (BaseInt bin = 0; bin < bins_data.bins.size(); ++bin) { - Cost max_load = model.bin_capacity() - item_weight; + Cost max_load = bin_capacity - item_weight; if (bins_data.loads[bin] <= max_load && (selected_bin == bins_data.bins.size() || bins_data.loads[bin] > bins_data.loads[selected_bin])) { @@ -232,9 +234,9 @@ void AddRandomizedBins(const BinPackingModel& model, BaseInt num_bins, // Generate bins all containing a specific item for (ElementIndex n : model.ItemRange()) { BaseInt unique_bin_num = scp_model.full_model().num_subsets(); - VLOG_EVERY_N_SEC(1, 5) - << "Generating bins: " << unique_bin_num << " / " << num_bins << " (" - << 100.0 * unique_bin_num / num_bins << "%)"; + VLOG_EVERY_N_SEC(1, 1) + << "[RGEN] Generating bins: " << unique_bin_num << " / " << num_bins + << " (" << 100.0 * unique_bin_num / num_bins << "%)"; if (scp_model.full_model().num_subsets() >= num_bins) { break; } @@ -251,7 +253,7 @@ void AddRandomizedBins(const BinPackingModel& model, BaseInt num_bins, bins_data.bins.push_back({n}); bins_data.loads.push_back(model.weights()[n]); } - BestFit(model, items, bins_data); + BestFit(model.weights(), model.bin_capacity(), items, bins_data); InsertBinsIntoModel(bins_data, scp_model); items.push_back(n); @@ -272,51 +274,86 @@ BinPackingSetCoverModel GenerateInitialBins(const BinPackingModel& model) { std::vector items(model.num_items()); absl::c_iota(items, ElementIndex(0)); - BestFit(model, items, bins_data); + absl::c_sort(items, [&](auto i1, auto i2) { + return model.weights()[i1] > model.weights()[i2]; + }); + BestFit(model.weights(), model.bin_capacity(), items, bins_data); InsertBinsIntoModel(bins_data, scp_model); - BaseInt solution_bin_num = bins_data.bins.size(); - VLOG(1) << "Best-fit solution: " << solution_bin_num << " bins"; - - // Largest first - if (!absl::c_is_sorted(model.weights(), std::greater<>())) { - bins_data.bins.clear(); - bins_data.loads.clear(); - absl::c_sort(items, [&](auto i1, auto i2) { - return model.weights()[i1] > model.weights()[i2]; - }); - BestFit(model, items, bins_data); - InsertBinsIntoModel(bins_data, scp_model); - } + VLOG(1) << "[BFIT] Largest first best-fit solution: " << bins_data.bins.size() + << " bins"; scp_model.CompleteModel(); return scp_model; } -void ExpKnap::Solve(const ElementCostVector& profits, - const ElementCostVector& weights, Cost capacity, - BaseInt bnb_nodes_limit) { +void ExpKnap::SaveBin() { + collected_bins_.back().clear(); + for (ElementIndex i : exceptions_) { + break_selection_[i] = !break_selection_[i]; + } + for (ElementIndex i : break_solution_) { + if (break_selection_[i]) { + collected_bins_.back().push_back(i); + } + } + for (ElementIndex i : exceptions_) { + if (break_selection_[i]) { + collected_bins_.back().push_back(i); + } + break_selection_[i] = !break_selection_[i]; + } + absl::c_sort(collected_bins_.back()); + DCHECK(absl::c_adjacent_find(collected_bins_.back()) == + collected_bins_.back().end()); +} +void ExpKnap::InitSolver(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit) { capacity_ = capacity; - best_delta_ = .0; - exceptions_.clear(); - maximal_exceptions_.clear(); - break_solution_.assign(profits.size(), false); items_.resize(profits.size()); - bnb_node_countdown_ = bnb_nodes_limit; - + collected_bins_.clear(); for (ElementIndex i; i < ElementIndex(items_.size()); ++i) { - items_[i] = {std::max(1e-6, profits[i]), weights[i], i}; + items_[i] = {profits[i], weights[i], i}; } - absl::c_sort(items_, [](Item i1, Item i2) { - return i1.profit / i1.weight > i2.profit / i2.weight; + return i1.profit / i1.weight < i2.profit / i2.weight; }); - Heuristic(items_); - maximal_exceptions_.push_back(exceptions_); + + bnb_node_countdown_ = bnb_nodes_limit; + best_delta_ = .0; exceptions_.clear(); - VLOG(5) << "[KPCG] Heuristic solution: cost " - << break_profit_sum_ + best_delta_; + break_selection_.assign(profits.size(), false); + break_solution_.clear(); +} + +void ExpKnap::FindGoodColumns(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit) { + InitSolver(profits, weights, capacity, bnb_nodes_limit); + Cost curr_best_cost = .0; + PartialBins more_bins; + std::vector remaining_items; - EleBranch(.0, break_weight_sum_ - capacity, break_it_ - 1, break_it_ + 1); + bnb_node_countdown_ = bnb_nodes_limit; + inserted_items_.assign(profits.size(), false); + do { + collected_bins_.emplace_back(); + Heuristic(); + VLOG(5) << "[KPCG] Heuristic solution: cost " + << break_profit_sum_ + best_delta_; + EleBranch(); + + for (ElementIndex i : collected_bins_.back()) { + inserted_items_[i] = true; + } + gtl::STLEraseAllFromSequenceIf( + &items_, [&](Item item) { return inserted_items_[item.index]; }); + } while (!items_.empty() && break_profit_sum_ + best_delta_ > 1.0); +} + +bool ExpKnap::EleBranch() { + exceptions_.clear(); + return EleBranch(.0, break_weight_sum_ - capacity_, break_it_ - 1, break_it_); } namespace { @@ -331,6 +368,8 @@ static Cost BoundCheck(Cost best_delta, Cost profit_delta, Cost overweight, // https://hjemmesider.diku.dk/~pisinger/expknap.c bool ExpKnap::EleBranch(Cost profit_delta, Cost overweigth, ItemIt out_item, ItemIt in_item) { + VLOG(6) << "[KPCG] EleBranch: profit_delta " << profit_delta << " overweigth " + << overweigth; if (bnb_node_countdown_-- <= 0) { return false; } @@ -345,7 +384,7 @@ bool ExpKnap::EleBranch(Cost profit_delta, Cost overweigth, ItemIt out_item, } bool maximal = true; - while (bnb_node_countdown_ > 0 && in_item < items_.end() && + while (bnb_node_countdown_ > 0 && in_item < items_.rend() && BoundCheck(best_delta_, profit_delta, overweigth, *in_item) >= 0) { exceptions_.push_back(in_item->index); Cost next_delta = profit_delta + in_item->profit; @@ -355,11 +394,11 @@ bool ExpKnap::EleBranch(Cost profit_delta, Cost overweigth, ItemIt out_item, } if (improved && maximal) { - maximal_exceptions_.push_back(exceptions_); + SaveBin(); } improved |= !maximal; } else { - while (bnb_node_countdown_ > 0 && out_item >= items_.begin() && + while (bnb_node_countdown_ > 0 && out_item >= items_.rbegin() && BoundCheck(best_delta_, profit_delta, overweigth, *out_item) >= 0) { exceptions_.push_back(out_item->index); Cost next_delta = profit_delta - out_item->profit; @@ -371,17 +410,18 @@ bool ExpKnap::EleBranch(Cost profit_delta, Cost overweigth, ItemIt out_item, return improved; } -void ExpKnap::Heuristic( - const util_intops::StrongVector& items) { +void ExpKnap::Heuristic() { + best_delta_ = break_profit_sum_ = break_weight_sum_ = .0; + break_it_ = items_.rbegin(); + break_selection_.assign(items_.size(), false); + break_solution_.clear(); exceptions_.clear(); - - break_profit_sum_ = break_weight_sum_ = .0; - break_it_ = items.begin(); - while (break_it_ < items.end() && + while (break_it_ < items_.rend() && break_it_->weight <= capacity_ - break_weight_sum_) { break_profit_sum_ += break_it_->profit; break_weight_sum_ += break_it_->weight; - break_solution_[break_it_->index] = true; + break_selection_[break_it_->index] = true; + break_solution_.push_back(break_it_->index); ++break_it_; } Cost residual = capacity_ - break_weight_sum_; @@ -391,16 +431,18 @@ void ExpKnap::Heuristic( Cost profit_delta_ub = residual * break_it_->profit / break_it_->weight; if (profit_delta_ub == .0) { + SaveBin(); return; } // Try filling the residual space with less efficient (maybe smaller) items best_delta_ = .0; - for (auto it = break_it_; it < items.end(); it++) { + for (auto it = break_it_; it < items_.rend(); it++) { if (it->weight <= residual && it->profit > best_delta_) { exceptions_ = {it->index}; best_delta_ = it->profit; if (best_delta_ >= profit_delta_ub) { + SaveBin(); return; } } @@ -408,57 +450,52 @@ void ExpKnap::Heuristic( // Try removing an item and adding the break item Cost min_weight = break_it_->weight - residual; - for (auto it = break_it_ - 1; it >= items.begin(); it--) { + for (auto it = break_it_ - 1; it >= items_.rbegin(); it--) { Cost profit_delta = break_it_->profit - it->profit; if (it->weight >= min_weight && profit_delta > best_delta_) { exceptions_ = {break_it_->index, it->index}; best_delta_ = profit_delta; if (best_delta_ >= profit_delta_ub) { + SaveBin(); return; } } } + SaveBin(); } bool BinPackingSetCoverModel::UpdateCore( Cost best_lower_bound, const ElementCostVector& best_multipliers, const scp::Solution& best_solution, bool force) { - if (--column_gen_countdown_ <= 0 && best_lower_bound != prev_lower_bound_) { - column_gen_countdown_ = column_gen_period_; + if (!base::IsTimeToUpdate(best_lower_bound, force)) { + return false; + } + + Cost full_lower_bound = base::UpdateMultipliers(best_multipliers); + if (scp::DivideIfGE0(std::abs(full_lower_bound - best_lower_bound), + best_lower_bound) < 0.01 && + best_lower_bound != prev_lower_bound_) { prev_lower_bound_ = best_lower_bound; - knapsack_solver_.Solve(best_multipliers, bpp_model_->weights(), - bpp_model_->bin_capacity(), - /*bnb_nodes_limit=*/10000); - const auto& exception_list = knapsack_solver_.maximal_exceptions(); - ElementBoolVector solution = knapsack_solver_.break_solution(); + knapsack_solver_.FindGoodColumns(best_multipliers, bpp_model_->weights(), + bpp_model_->bin_capacity(), + /*bnb_nodes_limit=*/1000); BaseInt num_added_bins = 0; - SparseColumn bin; - for (const std::vector& exception : exception_list) { - for (ElementIndex i : exception) { - solution[i] = !solution[i]; - } - - bin.clear(); - for (ElementIndex i : globals_.full_model.ElementRange()) { - if (solution[i]) { - bin.push_back(i); - } - } + for (const SparseColumn& bin : knapsack_solver_.collected_bins()) { num_added_bins += AddBin(bin) ? 1 : 0; - - for (ElementIndex i : exception) { - solution[i] = !solution[i]; - } } if (num_added_bins > 0) { VLOG(4) << "[KPCG] Added " << num_added_bins << " / " << globals_.full_model.num_subsets() << " bins"; } + base::SizeUpdate(); + // TODO(user): add incremental update only for the new columns just added + base::UpdateMultipliers(best_multipliers); } - scp::FullToCoreModel::SizeUpdate(); - scp::FullToCoreModel::FullToSubModelInvariantCheck(); - return scp::FullToCoreModel::UpdateCore(best_lower_bound, best_multipliers, - best_solution, force); + if (base::num_focus_subsets() < FixingFullModelView().num_focus_subsets()) { + ComputeAndSetFocus(best_lower_bound, best_solution); + return true; + } + return false; } } // namespace operations_research diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h index 0e627b6f5e0..17ce1417841 100644 --- a/ortools/algorithms/bin_packing.h +++ b/ortools/algorithms/bin_packing.h @@ -21,6 +21,7 @@ #include #include +#include #include #include "ortools/set_cover/base_types.h" @@ -61,21 +62,26 @@ class ExpKnap { Cost weight; // weight ElementIndex index; }; - using ItemIt = std::vector::const_iterator; + using ItemIt = std::vector::const_reverse_iterator; - void Solve(const ElementCostVector& profits, const ElementCostVector& weights, - Cost capacity, BaseInt bnb_nodes_limit); + void SaveBin(); + void FindGoodColumns(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit); + void InitSolver(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit); + bool EleBranch(); bool EleBranch(Cost profit_sum, Cost overweight, ItemIt out_item, ItemIt in_item); - void Heuristic(const util_intops::StrongVector& items); + void Heuristic(); - ElementBoolVector break_solution() const { return break_solution_; } - - std::vector> maximal_exceptions() const { - return maximal_exceptions_; + const std::vector& collected_bins() const { + return collected_bins_; } + Cost best_cost() const { return break_profit_sum_ + best_delta_; } private: Cost capacity_; // capacity @@ -85,12 +91,17 @@ class ExpKnap { Cost break_weight_sum_; Cost best_delta_; std::vector exceptions_; - std::vector> maximal_exceptions_; - ElementBoolVector break_solution_; + std::vector collected_bins_; + ElementBoolVector break_selection_; + ElementBoolVector inserted_items_; + SparseColumn break_solution_; BaseInt bnb_node_countdown_; + std::mt19937 rnd_; }; class BinPackingSetCoverModel : public scp::FullToCoreModel { + using base = scp::FullToCoreModel; + struct BinPackingModelGlobals { // Dirty hack to avoid invalidation of pointers/references // A pointer to this data structure is used to compute the hash of bins From 7cb1c36d2c8b095ae216333da0361b88ea13f55e Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Thu, 1 May 2025 18:31:52 +0200 Subject: [PATCH 35/35] Added small test to bin_packing example --- ortools/algorithms/samples/bin_packing_cft.cc | 87 +++++++++++++++++-- 1 file changed, 80 insertions(+), 7 deletions(-) diff --git a/ortools/algorithms/samples/bin_packing_cft.cc b/ortools/algorithms/samples/bin_packing_cft.cc index 2b297c6cf6b..0eafd474160 100644 --- a/ortools/algorithms/samples/bin_packing_cft.cc +++ b/ortools/algorithms/samples/bin_packing_cft.cc @@ -13,34 +13,107 @@ // limitations under the License. #include +#include #include #include #include "ortools/algorithms/bin_packing.h" #include "ortools/base/init_google.h" +#include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_cft.h" using namespace operations_research; ABSL_FLAG(std::string, instance, "", "BPP instance int RAIL format."); ABSL_FLAG(int, bins, 1000, "Number of bins to generate."); +template +std::string Stringify(const Iterable& col) { + std::string result; + for (auto i : col) { + absl::StrAppend(&result, " ", i); + } + return result; +} + +bool operator==(const SparseColumn& lhs, const std::vector& rhs) { + if (lhs.size() != rhs.size()) return false; + auto lit = lhs.begin(); + auto rit = rhs.begin(); + while (lit != lhs.end() && rit != rhs.end()) { + if (static_cast(*lit) != *rit) { + return false; + } + ++lit; + ++rit; + } + return true; +} + +void RunTest(const ElementCostVector& weights, const ElementCostVector& profits, + const std::vector& expected) { + ExpKnap knap_solver; + + for (ElementIndex i; i < ElementIndex(weights.size()); ++i) { + std::cout << "Item " << i << " -- profit: " << profits[i] + << " weight: " << weights[i] + << " efficiency: " << profits[i] / weights[i] << "\n"; + } + + knap_solver.InitSolver(profits, weights, 6, 100000000); + knap_solver.Heuristic(); + std::cout << "Heur solution cost " << knap_solver.best_cost() << " -- " + << Stringify(knap_solver.collected_bins()[0]) << "\n"; + + knap_solver.EleBranch(); + std::cout << "B&b solution cost " << knap_solver.best_cost() << " -- " + << Stringify(knap_solver.collected_bins()[0]) << "\n"; + + const auto& result = knap_solver.collected_bins()[0]; + if (!(result == expected)) { + std::cout << "Error: expected " << Stringify(expected) << " but got " + << Stringify(result) << "\n"; + } + std::cout << std::endl; +} + +void KnapsackTest() { + std::cout << "Testing knapsack\n"; + ExpKnap knap_solver; + ElementCostVector ws = {1, 2, 3, 4, 5}; + RunTest(ws, {10, 20, 30, 40, 51}, {0, 4}); + RunTest(ws, {10, 20, 30, 41, 50}, {1, 3}); + RunTest(ws, {10, 20, 31, 40, 50}, {0, 1, 2}); + RunTest(ws, {10, 21, 30, 41, 50}, {1, 3}); + RunTest(ws, {11, 21, 30, 40, 50}, {0, 1, 2}); + RunTest(ws, {11, 20, 31, 40, 50}, {0, 1, 2}); + RunTest(ws, {11, 20, 30, 41, 50}, {0, 4}); + RunTest(ws, {11, 20, 30, 40, 51}, {0, 4}); + RunTest(ws, {11, 21, 31, 40, 50}, {0, 1, 2}); + RunTest({4.1, 2, 2, 2}, {8.5, 3, 3, 3}, {1, 2, 3}); +} + int main(int argc, char** argv) { InitGoogle(argv[0], &argc, &argv, true); + // KnapsackTest(); + // return 0; + BinPackingModel model = ReadBpp(absl::GetFlag(FLAGS_instance)); // Quick run with a minimal set of bins BinPackingSetCoverModel scp_model = GenerateInitialBins(model); scp::PrimalDualState best_result = scp::RunCftHeuristic(scp_model); - // Run the CFT again with more bins to get a better solution - std::mt19937 rnd(0); - AddRandomizedBins(model, absl::GetFlag(FLAGS_bins), scp_model, rnd); - scp::PrimalDualState result = - scp::RunCftHeuristic(scp_model, best_result.solution); - if (result.solution.cost() < best_result.solution.cost()) { - best_result = result; + if (absl::GetFlag(FLAGS_bins) > 0) { + // Run the CFT again with more bins to get a better solution + std::mt19937 rnd(0); + AddRandomizedBins(model, absl::GetFlag(FLAGS_bins), scp_model, rnd); + scp::PrimalDualState result = + scp::RunCftHeuristic(scp_model, best_result.solution); + if (result.solution.cost() < best_result.solution.cost()) { + best_result = result; + } } auto [solution, dual] = best_result;