diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc new file mode 100644 index 00000000000..71628ed1b30 --- /dev/null +++ b/ortools/set_cover/set_cover_cft.cc @@ -0,0 +1,1038 @@ +// 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/set_cover/set_cover_cft.h" + +#include +#include +#include +#include +#include + +#include +#include + +#include "ortools/base/status_macros.h" +#include "ortools/base/stl_util.h" +#include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_views.h" + +namespace operations_research::scp { + +//////////////////////////////////////////////////////////////////////// +////////////////////////// COMMON DEFINITIONS ////////////////////////// +//////////////////////////////////////////////////////////////////////// +namespace { +class CoverCounters { + public: + CoverCounters(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())); + return cov_counters[i]; + } + + template + BaseInt Cover(const IterableT& subset) { + BaseInt covered = 0; + for (ElementIndex i : subset) { + covered += cov_counters[i] == 0 ? 1ULL : 0ULL; + cov_counters[i]++; + } + return covered; + } + + template + BaseInt Uncover(const IterableT& subset) { + BaseInt uncovered = 0; + for (ElementIndex i : subset) { + --cov_counters[i]; + uncovered += cov_counters[i] == 0 ? 1ULL : 0ULL; + } + return uncovered; + } + + // Check if all the elements of a subset are already covered. + template + bool IsRedundantCover(IterableT const& subset) const { + return absl::c_all_of(subset, + [&](ElementIndex 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; }); + } + + private: + ElementToIntVector cov_counters; +}; +} // namespace + +Solution::Solution(const SubModel& model, + const std::vector& core_subsets) + : cost_(), subsets_() { + subsets_.reserve(core_subsets.size() + model.fixed_columns().size()); + cost_ = model.fixed_cost(); + for (FullSubsetIndex full_j : model.fixed_columns()) { + subsets_.push_back(full_j); + } + for (SubsetIndex core_j : core_subsets) { + FullSubsetIndex full_j = model.MapCoreToFullSubsetIndex(core_j); + AddSubset(full_j, model.subset_costs()[core_j]); + } +} + +/////////////////////////////////////////////////////////////////////// +///////////////////////////// SUBGRADIENT ///////////////////////////// +/////////////////////////////////////////////////////////////////////// + +BoundCBs::BoundCBs(const SubModel& model) + : squared_norm_(static_cast(model.num_elements())), + direction_(ElementCostVector(model.num_elements(), .0)), + 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_period_(300), // Arbitrary from [1] + 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] + step_size_update_period_(20) // Arbitrary from [1] +{} + +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) { + return true; + } + if (--exit_test_countdown_ > 0) { + return false; + } + exit_test_countdown_ = exit_test_period_; + Cost best_lower_bound = context.best_dual_state.lower_bound(); + Cost abs_improvement = best_lower_bound - prev_best_lb_; + Cost rel_improvement = DivideIfGE0(abs_improvement, best_lower_bound); + prev_best_lb_ = best_lower_bound; + return abs_improvement < 1.0 && rel_improvement < .001; +} + +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_; + } + + squared_norm_ = .0; + // Stabilize current direction and compute squared norm + for (ElementIndex i : context.model.ElementRange()) { + direction_[i] = stabilization_coeff * direction_[i] + + (1.0 - stabilization_coeff) * prev_direction_[i]; + if ((context.current_dual_state.multipliers()[i] <= .0) && + (direction_[i] < .0)) { + direction_[i] = 0; + } + squared_norm_ += direction_[i] * direction_[i]; + prev_direction_[i] = direction_[i]; + } + + 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]; + } +} + +void BoundCBs::UpdateStepSize(Cost 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); + + if (--step_size_update_countdown_ <= 0) { + step_size_update_countdown_ = step_size_update_period_; + + 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; + + } else if (gap > .01) { // Arbitray from [1] + step_size_ /= 2.0; // Arbitray from [1] + } + last_min_lb_seen_ = std::numeric_limits::max(); + last_max_lb_seen_ = .0; + // Not described in the paper, but in rare cases the subgradient diverges + step_size_ = std::clamp(step_size_, 1e-6, 10.0); // Arbitrary from c4v4 + } +} + +void BoundCBs::MakeMinimalCoverageSubgradient(const SubgradientContext& context, + ElementCostVector& subgradient) { + lagrangian_solution_.clear(); + const auto& reduced_costs = context.current_dual_state.reduced_costs(); + for (SubsetIndex j : context.model.SubsetRange()) { + if (reduced_costs[j] < .0) { + lagrangian_solution_.push_back(j); + } + } + + absl::c_sort(lagrangian_solution_, [&](SubsetIndex j1, SubsetIndex j2) { + return reduced_costs[j1] > reduced_costs[j2]; + }); + + const auto& cols = context.model.columns(); + for (SubsetIndex j : lagrangian_solution_) { + if (absl::c_all_of(cols[j], [&](auto i) { return subgradient[i] < 0; })) { + absl::c_for_each(cols[j], [&](auto i) { subgradient[i] += 1.0; }); + } + } +} + +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; + return true; + } + return false; +} + +void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, + PrimalDualState& best_state) { + DCHECK(ValidateSubModel(model)); + + ElementCostVector subgradient = ElementCostVector(model.num_elements(), .0); + DualState dual_state = best_state.dual_state; + 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_solution = best_state.solution, + .subgradient = subgradient}; + size_t iter = 0; + while (!cbs.ExitCondition(context)) { + // Compute subgradient (O(nnz)) + subgradient.assign(model.num_elements(), 1.0); + for (SubsetIndex j : model.SubsetRange()) { + if (dual_state.reduced_costs()[j] < .0) { + for (ElementIndex i : model.columns()[j]) { + subgradient[i] -= 1.0; + } + } + } + + cbs.ComputeMultipliersDelta(context, multipliers_delta); + dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { + i_mult = std::max(.0, i_mult + multipliers_delta[i]); + }); + if (dual_state.lower_bound() > best_state.dual_state.lower_bound()) { + best_state.dual_state = dual_state; + } + + cbs.RunHeuristic(context, solution); + if (!solution.subsets().empty() && + solution.cost() < best_state.solution.cost()) { + best_state.solution = solution; + } + + if (++iter % 50 == 0) + std::cout << "Subgradient " << iter << " -- Bounds: Lower " + << dual_state.lower_bound() << ", best " + << best_state.dual_state.lower_bound() << " - Upper " + << best_state.solution.cost() - model.fixed_cost() + << ", global " << best_state.solution.cost() << "\n"; + + if (cbs.UpdateCoreModel(model, best_state)) { + std::cout << "Core model has been updated.\n"; + dual_state = best_state.dual_state; + } + } + 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"; +} + +//////////////////////////////////////////////////////////////////////// +/////////////////////// MULTIPLIERS BASED GREEDY /////////////////////// +//////////////////////////////////////////////////////////////////////// +namespace { +struct Score { + Cost score; + SubsetIndex idx; +}; + +class GreedyScores { + public: + static constexpr BaseInt removed_idx = -1; + static constexpr Cost max_score = std::numeric_limits::max(); + + GreedyScores(const SubModel& model, const DualState& dual_state) + : bad_size_(), + worst_good_score_(std::numeric_limits::lowest()), + scores_(), + reduced_costs_(dual_state.reduced_costs()), + covering_counts_(model.num_subsets()), + score_map_(model.num_subsets()) { + BaseInt s = 0; + for (SubsetIndex j : model.SubsetRange()) { + DCHECK(model.column_size(j) > 0); + covering_counts_[j] = model.column_size(j); + Cost j_score = ComputeScore(reduced_costs_[j], covering_counts_[j]); + scores_.push_back({j_score, j}); + score_map_[j] = s++; + DCHECK(std::isfinite(reduced_costs_[j])); + DCHECK(std::isfinite(j_score)); + } + bad_size_ = scores_.size(); + } + + SubsetIndex FindMinScoreColumn(const SubModel& model, + const DualState& dual_state) { + // Check if the bad/good partition should be updated + if (bad_size_ == scores_.size()) { + if (bad_size_ > model.num_focus_elements()) { + bad_size_ = bad_size_ - model.num_focus_elements(); + absl::c_nth_element(scores_, scores_.begin() + bad_size_, + [](Score a, Score b) { return a.score > b.score; }); + worst_good_score_ = scores_[bad_size_].score; + for (BaseInt s = 0; s < scores_.size(); ++s) { + score_map_[scores_[s].idx] = s; + } + } else { + bad_size_ = 0; + worst_good_score_ = max_score; + } + DCHECK(bad_size_ > 0 || worst_good_score_ == max_score); + } + + Score min_score = + *std::min_element(scores_.begin() + bad_size_, scores_.end(), + [](Score a, Score b) { return a.score < b.score; }); + SubsetIndex j_star = min_score.idx; + DCHECK_LT(min_score.score, max_score); + return j_star; + } + + // For each row in the given set, if `cond` returns true, the row is + // considered newly covered. The function then iterates over the columns of + // that row, updating the scores of the columns accordingly. + template + BaseInt UpdateColumnsScoreOfRowsIf(const RowT& rows, + const ElementCostVector& multipliers, + const ElementSpanT& row_idxs, CondT cond) { + BaseInt processed_rows_count = 0; + for (ElementIndex i : row_idxs) { + if (!cond(i)) { + continue; + } + + ++processed_rows_count; + for (SubsetIndex j : rows[i]) { + covering_counts_[j] -= 1; + reduced_costs_[j] += multipliers[i]; + + BaseInt s = score_map_[j]; + DCHECK_NE(s, removed_idx) << "Column is not in the score map"; + scores_[s].score = ComputeScore(reduced_costs_[j], covering_counts_[j]); + + if (covering_counts_[j] == 0) { + // Column is redundant: its score can be removed + if (s < bad_size_) { + // Column is bad: promote to good partition before removal + SwapScores(s, bad_size_ - 1); + s = --bad_size_; + } + SwapScores(s, scores_.size() - 1); + scores_.pop_back(); + } else if (s >= bad_size_ && scores_[s].score > worst_good_score_) { + // Column not good anymore: move it into bad partition + SwapScores(s, bad_size_++); + } + } + } + return processed_rows_count; + } + + private: + void SwapScores(BaseInt s1, BaseInt s2) { + SubsetIndex j1 = scores_[s1].idx, j2 = scores_[s2].idx; + std::swap(scores_[s1], scores_[s2]); + std::swap(score_map_[j1], score_map_[j2]); + } + + // Score computed as described in [1] + static Cost ComputeScore(Cost adjusted_reduced_cost, + BaseInt num_rows_covered) { + DCHECK(std::isfinite(adjusted_reduced_cost)) << "Gamma is not finite"; + return num_rows_covered == 0 + ? max_score + : (adjusted_reduced_cost > .0 + ? adjusted_reduced_cost / num_rows_covered + : adjusted_reduced_cost * num_rows_covered); + } + + private: + // scores_ is partitioned into bad-scores / good-scores + BaseInt bad_size_; + + // sentinel level to trigger a partition update of the scores + Cost worst_good_score_; + + // column scores kept updated + std::vector scores_; + + // reduced costs adjusted to currently uncovered rows (size=n) + SubsetCostVector reduced_costs_; + + // number of uncovered rows covered by each column (size=n) + SubsetToIntVector covering_counts_; + + // position of each column score into the scores_ + SubsetToIntVector score_map_; +}; + +// Stores the redundancy set and related information +class RedundancyRemover { + public: + RedundancyRemover(const SubModel& model, CoverCounters& total_coverage) + : redund_set_(), + total_coverage_(total_coverage), + partial_coverage_(model.num_elements()), + partial_cost_(.0), + partial_size_(0), + partial_cov_count_(0), + cols_to_remove_() {} + + Cost TryRemoveRedundantCols(const SubModel& model, Cost cost_cutoff, + std::vector& sol_subsets) { + for (SubsetIndex j : sol_subsets) { + if (total_coverage_.IsRedundantUncover(model.columns()[j])) + redund_set_.push_back({model.subset_costs()[j], j}); + else { + ++partial_size_; + partial_cost_ += model.subset_costs()[j]; + partial_cov_count_ += partial_coverage_.Cover(model.columns()[j]); + } + if (partial_cost_ >= cost_cutoff) { + return partial_cost_; + } + } + if (redund_set_.empty()) { + return partial_cost_; + } + absl::c_sort(redund_set_, + [](Score a, Score b) { return a.score < b.score; }); + + if (partial_cov_count_ < model.num_focus_elements()) { + // Complete partial solution heuristically + HeuristicRedundancyRemoval(model, cost_cutoff); + } else { + // All redundant columns can be removed + for (Score redund_col : redund_set_) { + cols_to_remove_.push_back(redund_col.idx); + } + } + + // Note: In [1], an enumeration to selected the best redundant columns to + // remove is performed when the number of redundant columns is <= 10. + // However, based on experiments with github.com/c4v4/cft/, it appears + // that this enumeration does not provide significant benefits to justify + // the added complexity. + + if (partial_cost_ < cost_cutoff) { + gtl::STLEraseAllFromSequenceIf(&sol_subsets, [&](auto j) { + return absl::c_any_of(cols_to_remove_, + [j](auto j_r) { return j_r == j; }); + }); + } + return partial_cost_; + } + + private: + // Remove redundant columns from the redundancy set using a heuristic + void HeuristicRedundancyRemoval(const SubModel& model, Cost cost_cutoff) { + while (!redund_set_.empty()) { + if (partial_cov_count_ == model.num_focus_elements()) { + return; + } + SubsetIndex j = redund_set_.back().idx; + const auto& j_col = model.columns()[j]; + redund_set_.pop_back(); + + if (total_coverage_.IsRedundantUncover(j_col)) { + total_coverage_.Uncover(j_col); + cols_to_remove_.push_back(j); + } else { + partial_cost_ += model.subset_costs()[j]; + partial_cov_count_ += partial_coverage_.Cover(j_col); + } + } + } + + private: + // redundant columns + their cost + std::vector redund_set_; + + // row-cov if all the remaining columns are selected + CoverCounters total_coverage_; + + // row-cov if we selected the current column + CoverCounters partial_coverage_; + + // current partial solution cost + Cost partial_cost_; + + // current partial solution size + BaseInt partial_size_; + + // number of covered rows + Cost partial_cov_count_; + + // list of columns to remove + std::vector cols_to_remove_; +}; + +} // namespace + +Solution RunMultiplierBasedGreedy(const SubModel& model, + const DualState& dual_state, + Cost cost_cutoff) { + std::vector sol_subsets; + CoverGreedly(model, dual_state, cost_cutoff, + std::numeric_limits::max(), sol_subsets); + return Solution(model, sol_subsets); +} + +Cost CoverGreedly(const SubModel& model, const DualState& dual_state, + Cost cost_cutoff, BaseInt stop_size, + std::vector& sol_subsets) { + Cost sol_cost = .0; + for (SubsetIndex j : sol_subsets) { + sol_cost += model.subset_costs()[j]; + } + if (sol_cost >= cost_cutoff) { + sol_subsets.clear(); + return std::numeric_limits::max(); + } + if (sol_subsets.size() >= stop_size) { + // Solution already has required size -> early exit + return sol_cost; + } + + // Process input solution (if not empty) + BaseInt num_rows_to_cover = model.num_focus_elements(); + CoverCounters covered_rows(model.num_elements()); + for (SubsetIndex j : sol_subsets) { + num_rows_to_cover -= covered_rows.Cover(model.columns()[j]); + if (num_rows_to_cover == 0) { + return sol_cost; + } + } + + // Initialize column scores taking into account rows already covered + GreedyScores scores(model, dual_state); // TODO(?): cache it! + if (!sol_subsets.empty()) { + scores.UpdateColumnsScoreOfRowsIf( + model.rows(), dual_state.multipliers(), model.ElementRange(), + [&](auto i) { return covered_rows[i] > 0; }); + } + + // Fill up partial solution + while (num_rows_to_cover > 0 && sol_subsets.size() < stop_size) { + SubsetIndex j_star = scores.FindMinScoreColumn(model, dual_state); + const auto& column = model.columns()[j_star]; + num_rows_to_cover -= scores.UpdateColumnsScoreOfRowsIf( + model.rows(), dual_state.multipliers(), column, + [&](auto i) { return covered_rows[i] == 0; }); + sol_subsets.push_back(j_star); + covered_rows.Cover(column); + } + + // Either remove redundant columns or discard solution + RedundancyRemover remover(model, covered_rows); // TODO(?): cache it! + return remover.TryRemoveRedundantCols(model, cost_cutoff, sol_subsets); +} + +/////////////////////////////////////////////////////////////////////// +//////////////////////// THREE PHASE ALGORITHM //////////////////////// +/////////////////////////////////////////////////////////////////////// +namespace { + +DualState MakeTentativeDualState(const SubModel& model) { + DualState tentative_dual_state(model); + tentative_dual_state.DualUpdate( + model, [&](ElementIndex i, Cost& i_multiplier) { + i_multiplier = std::numeric_limits::max(); + for (SubsetIndex j : model.rows()[i]) { + Cost candidate = model.subset_costs()[j] / model.column_size(j); + i_multiplier = std::min(i_multiplier, candidate); + } + }); + return tentative_dual_state; +} + +void FixBestColumns(SubModel& model, PrimalDualState& state) { + auto& [best_sol, dual_state] = state; + + // This approach is not the most efficient but prioritizes clarity and the + // current abstraction system. We save the current core multipliers, mapped to + // the full model's element indices. By organizing the multipliers using the + // full model indices, we can easily map them to the new core model indices + // after fixing columns. Note: This mapping isn't strictly necessary because + // fixing columns effectively removes rows, and the remaining multipliers + // naturally shift to earlier positions. A simple iteration would suffice to + // discard multipliers for rows no longer in the new core model. + FullElementCostVector full_multipliers(model.num_elements(), .0); + for (ElementIndex core_i : model.ElementRange()) { + FullElementIndex full_i = model.MapCoreToFullElementIndex(core_i); + full_multipliers[full_i] = dual_state.multipliers()[core_i]; + } + + std::vector cols_to_fix; + CoverCounters row_coverage(model.num_elements()); + for (SubsetIndex j : model.SubsetRange()) { + if (dual_state.reduced_costs()[j] < -0.001) { + cols_to_fix.push_back(j); + row_coverage.Cover(model.columns()[j]); + } + } + + // Remove columns that overlap between each other + gtl::STLEraseAllFromSequenceIf(&cols_to_fix, [&](SubsetIndex j) { + return absl::c_any_of(model.columns()[j], + [&](ElementIndex i) { return row_coverage[i] > 1; }); + }); + + // Ensure at least a minimum number of columns are fixed + BaseInt fix_at_least = + cols_to_fix.size() + std::max(1, model.num_elements() / 200); + CoverGreedly(model, state.dual_state, std::numeric_limits::max(), + fix_at_least, cols_to_fix); + + // Fix columns and update the model + Cost fixed_cost_delta = model.FixColumns(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'; + + // Update multipliers for the reduced model + dual_state.DualUpdate(model, [&](ElementIndex core_i, Cost& multiplier) { + // Note: if SubModelView is used as CoreModel, then this mapping is always + // the identity mapping and can be removed + multiplier = full_multipliers[model.MapCoreToFullElementIndex(core_i)]; + }); +} + +void RandomizeDualState(const SubModel& model, DualState& dual_state, + std::mt19937& rnd) { + dual_state.DualUpdate(model, [&](ElementIndex, Cost& i_multiplier) { + i_multiplier *= absl::Uniform(rnd, 0.9, 1.1); + }); +} +} // namespace + +void HeuristicCBs::RunHeuristic(const SubgradientContext& context, + Solution& solution) { + solution = RunMultiplierBasedGreedy( + context.model, context.current_dual_state, + context.best_solution.cost() - context.model.fixed_cost()); +} + +void HeuristicCBs::ComputeMultipliersDelta(const SubgradientContext& context, + ElementCostVector& delta_mults) { + Cost squared_norm = .0; + for (ElementIndex i : context.model.ElementRange()) { + squared_norm += context.subgradient[i] * context.subgradient[i]; + } + + Cost lower_bound = context.current_dual_state.lower_bound(); + Cost upper_bound = context.best_solution.cost() - context.model.fixed_cost(); + DCHECK_GE(upper_bound, lower_bound); + 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 * context.subgradient[i]; + } +} + +PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { + DCHECK(ValidateSubModel(model)); + + PrimalDualState best_state = {.solution = init_solution, + .dual_state = MakeTentativeDualState(model)}; + if (best_state.solution.Empty()) { + 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"; + + 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) { + ++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'; + + // Phase 1: refine the current dual_state and model + BoundCBs dual_bound_cbs(model); + std::cout << "\nSubgradient Phase:\n"; + SubgradientOptimization(model, dual_bound_cbs, curr_state); + if (iter_count == 1) { + best_state.dual_state = curr_state.dual_state; + } + // Phase 2: search for good solutions + HeuristicCBs heuristic_cbs; + heuristic_cbs.set_step_size(dual_bound_cbs.step_size()); + std::cout << "\nHeuristic Phase:\n"; + SubgradientOptimization(model, heuristic_cbs, curr_state); + if (curr_state.solution.cost() < best_state.solution.cost()) { + best_state.solution = curr_state.solution; + } + + std::cout << "\n3Phase Bounds: Lower " + << best_state.dual_state.lower_bound() << ", Upper " + << best_state.solution.cost() << '\n'; + + // Phase 3: Fix the best columns (diving) + FixBestColumns(model, curr_state); + RandomizeDualState(model, curr_state.dual_state, rnd); + } + + std::cout << "\n\n3Phase End\n"; + std::cout << "Final Bounds: Lower " << best_state.dual_state.lower_bound() + << ", Upper " << best_state.solution.cost() << '\n'; + + return best_state; +} + +/////////////////////////////////////////////////////////////////////// +//////////////////////// FULL TO CORE PRICING ///////////////////////// +/////////////////////////////////////////////////////////////////////// + +namespace { +std::vector ComputeTentativeFocus(StrongModelView full_model) { + FullSubsetBoolVector selected(full_model.num_subsets(), false); + std::vector columns_focus; + columns_focus.reserve(full_model.num_elements() * kMinCov); + + // Select the first min_row_coverage columns for each row + for (const auto& row : full_model.rows()) { + BaseInt countdown = kMinCov; + for (FullSubsetIndex j : row) { + if (--countdown <= 0) { + break; + } + if (!selected[j]) { + selected[j] = true; + columns_focus.push_back(j); + } + } + } + // NOTE: unnecessary, but it keeps equivalence between SubModelView/SubModel + absl::c_sort(columns_focus); + return columns_focus; +} + +void SelecteMinRedCostColumns(FilterModelView full_model, + const SubsetCostVector& reduced_costs, + std::vector& new_core_columns, + FullSubsetBoolVector& selected) { + DCHECK_EQ(reduced_costs.size(), full_model.num_subsets()); + DCHECK_EQ(selected.size(), full_model.num_subsets()); + + std::vector candidates; + for (SubsetIndex j : full_model.SubsetRange()) + if (reduced_costs[j] < 0.1) { + candidates.push_back(j); + } + + BaseInt max_size = 5 * full_model.num_focus_elements(); + if (candidates.size() > max_size) { + absl::c_nth_element(candidates, candidates.begin() + max_size - 1, + [&](SubsetIndex j1, SubsetIndex j2) { + return reduced_costs[j1] < reduced_costs[j2]; + }); + candidates.resize(max_size); + } + for (SubsetIndex j : candidates) { + FullSubsetIndex j_full = static_cast(j); + if (!selected[j_full]) { + selected[j_full] = true; + new_core_columns.push_back(j_full); + } + } +} + +static void SelectMinRedCostByRow(FilterModelView full_model, + const SubsetCostVector& reduced_costs, + std::vector& columns_map, + FullSubsetBoolVector& selected) { + DCHECK_EQ(reduced_costs.size(), full_model.num_subsets()); + DCHECK_EQ(selected.size(), full_model.num_subsets()); + + for (const auto& row : full_model.rows()) { + // Collect best `kMinCov` columns covering row `i` + SubsetIndex best_cols[kMinCov]; + BaseInt best_size = 0; + for (SubsetIndex j : row) { + if (best_size < kMinCov) { + best_cols[best_size++] = j; + continue; + } + if (reduced_costs[j] < reduced_costs[best_cols[kMinCov - 1]]) { + BaseInt n = kMinCov - 1; + while (n > 0 && reduced_costs[j] < reduced_costs[best_cols[n - 1]]) { + best_cols[n] = best_cols[n - 1]; + --n; + } + best_cols[n] = j; + } + } + + DCHECK(best_size > 0); + for (BaseInt s = 0; s < best_size; ++s) { + FullSubsetIndex j = static_cast(best_cols[s]); + if (!selected[j]) { + selected[j] = true; + columns_map.push_back(j); + } + } + } +} +} // namespace + +FullToCoreModel::FullToCoreModel(const Model* full_model) + : SubModel(full_model, ComputeTentativeFocus(StrongModelView(full_model))), + full_model_(full_model), + is_focus_col_(full_model->num_subsets(), true), + is_focus_row_(full_model->num_elements(), true), + num_subsets_(full_model->num_subsets()), + num_elements_(full_model->num_elements()), + full_dual_state_(*full_model), + best_dual_state_(full_dual_state_) { + ResetPricingPeriod(); + DCHECK(ValidateSubModel(*this)); + DCHECK(FullToSubModelInvariantCheck()); +} + +void FullToCoreModel::ResetPricingPeriod() { + update_countdown_ = 10; + update_period_ = 10; + update_max_period_ = std::min(1000, full_model_->num_elements() / 3); +} + +Cost FullToCoreModel::FixColumns( + 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); + 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(); + Cost fixed_cost = base::FixColumns(columns_to_fix); + DCHECK(FullToSubModelInvariantCheck()); + return fixed_cost; +} + +bool FullToCoreModel::UpdateCore(PrimalDualState& core_state) { + if (--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'; + + 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()) { + if (IsFocusCol(full_j)) { + new_core_columns.push_back(full_j); + selected_columns[full_j] = true; + } + } + + 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); + + // NOTE: unnecessary, but it keeps equivalence between SubModelView/SubModel + absl::c_sort(new_core_columns); + SetFocus(new_core_columns); + core_state.dual_state.DualUpdate(*this, [](ElementIndex i, Cost& i_mult) { + // multipliers didn't cange, but reduced cost must be recomputed + }); + + 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()); + if (ratio <= 1e-6) { + update_period_ = std::min(update_max_period_, 10 * update_period_); + } else if (ratio <= 0.02) { + update_period_ = std::min(update_max_period_, 5 * update_period_); + } else if (ratio <= 0.2) { + update_period_ = std::min(update_max_period_, 2 * update_period_); + } else { + update_period_ = 10; + } + update_countdown_ = update_period_; +} + +void FullToCoreModel::UpdateDualState(const DualState& core_dual_state, + 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]; + }); + + // Here, we simply check if any columns have been fixed, and only update the + // best dual state when no fixing is in place. + // + // Mapping a "local" dual state to a global one is possible. This would + // involve keeping the multipliers for non-focused elements fixed, updating + // the multipliers for focused elements, and then computing the dual state for + // the entire model. However, this approach is not implemented here. Such a + // strategy is unlikely to improve the best dual state unless the focus is + // *always* limited to a small subset of elements (and therefore the LB sucks + // and it is easy to improve) and the CFT is applied exclusively within that + // narrow scope, but this falls outside the current scope of this project. + if (fixed_columns().empty() && + full_dual_state_.lower_bound() > best_dual_state_.lower_bound()) { + best_dual_state_ = full_dual_state_; + } +} + +bool FullToCoreModel::FullToSubModelInvariantCheck() { + const SubModel& sub_model = *this; + StrongModelView typed_full_model = StrongTypedFullModelView(); + + if (typed_full_model.num_subsets() < sub_model.num_subsets()) { + std::cerr << absl::StrCat("SubModelView has ", sub_model.num_subsets(), + " subsets, but the full model has ", + typed_full_model.num_subsets(), " subsets.\n"); + return false; + } + if (typed_full_model.num_elements() != sub_model.num_elements()) { + std::cerr << absl::StrCat("SubModelView has ", sub_model.num_elements(), + " elements, but the full model has ", + typed_full_model.num_elements(), " elements.\n"); + return false; + } + for (SubsetIndex core_j : sub_model.SubsetRange()) { + FullSubsetIndex full_j = sub_model.MapCoreToFullSubsetIndex(core_j); + if (!is_focus_col_[static_cast(full_j)]) { + std::cerr << absl::StrCat("Subset ", core_j, + " in sub-model but its mapped subset ", full_j, + " not found 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)]) { + std::cerr << absl::StrCat("Element ", core_i, + " in sub-model but its mapped element ", full_i, + " not found in full model view.\n"); + return false; + } + } + for (FullElementIndex full_i : typed_full_model.ElementRange()) { + if (!IsFocusRow(full_i)) { + continue; + } + ElementIndex core_i = sub_model.MapFullToCoreElementIndex(full_i); + if (core_i < ElementIndex() || + ElementIndex(sub_model.num_elements()) < core_i) { + std::cerr << absl::StrCat( + "Element ", full_i, + " in full model view but has no mapped element in sub-model.\n"); + return false; + } + } + return true; +} + +} // namespace operations_research::scp \ No newline at end of file diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h new file mode 100644 index 00000000000..5cdc99123eb --- /dev/null +++ b/ortools/set_cover/set_cover_cft.h @@ -0,0 +1,407 @@ +// 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_SET_COVER_SET_COVER_CFT_H +#define OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H + +#include +#include +#include + +#include + +#include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_submodel.h" +#include "ortools/set_cover/set_cover_views.h" + +namespace operations_research::scp { + +// Implementation of: +// Caprara, Alberto, Matteo Fischetti, and Paolo Toth. 1999. “A Heuristic +// Method for the Set Covering Problem.” Operations Research 47 (5): 730–43. +// https://www.jstor.org/stable/223097 +// +// Hereafter referred to as CFT. +// +// SUMMARY +// The CFT algorithm is a heuristic approach to the set-covering problem. At its +// core, it combines a primal greedy heuristic with dual information obtained +// from the optimization of the Lagrangian relaxation of the problem. +// (Note: columns/subsets and rows/elements will be used interchangeably.) +// +// STRUCTURE +// The core of the algorithm is the 3Phase: +// 1. Subgradient optimization of the Lagrangian relaxation. +// 2. A primal greedy heuristic guided by the dual information. +// 3. Fixing some of the "best" columns (in terms of reduced costs) into the +// solution (diving). +// + Repeat until an exit criterion is met. +// +// The paper also considers an optional outer loop, which invokes the 3Phase +// process and then fixes some columns from the current best solution. This +// introduces two levels of diving: the outer loop fixes "primal-good" columns +// (based on the best solution), while the inner loop fixes "dual-good" columns +// (based on reduced costs). +// NOTE: The outer loop is not implemented in this version (yet - April 2025). +// +// Key characteristics of the algorithm: +// +// - The CFT algorithm is tailored for instances where the number of columns is +// significantly larger than the number of rows. +// +// - To improve efficiency, a core model approach is used. This involves +// selecting a small subset of columns based on their reduced costs, thereby +// substantially reducing the problem size handled in the internal steps. +// +// - Due to the use of the core model and column fixing, the algorithm rarely +// considers the entire problem. Instead, it operates on a small "window" of +// the problem. Efficiently managing this small window is a central aspect of +// any CFT implementation. +// +// - The core model scheme also enables an alternative implementation where the +// algorithm starts with a small model and progressively adds columns through +// a column-generation procedure. While this column generation procedure is +// problem-dependent and cannot be implemented here, the architecture of this +// implementation is designed to be extensible, allowing for such a procedure +// to be added in the future. +// + +//////////////////////////////////////////////////////////////////////// +////////////////////////// COMMON DEFINITIONS ////////////////////////// +//////////////////////////////////////////////////////////////////////// + +// Small class to store the solution of a sub-model. It contains the cost and +// the subset list. +class Solution { + public: + Solution() = default; + Solution(const SubModel& model, const std::vector& core_subsets); + + double cost() const { return cost_; } + const std::vector& subsets() const { return subsets_; } + void AddSubset(FullSubsetIndex subset, Cost cost) { + subsets_.push_back(subset); + cost_ += cost; + } + bool Empty() const { return subsets_.empty(); } + void Clear() { + cost_ = 0.0; + subsets_.clear(); + } + + private: + Cost cost_; + std::vector subsets_; +}; + +// In the narrow scope of the CFT subgradient, there are often divisions +// between non-negative quantities (e.g., to compute a relative gap). In these +// specific cases, the denominator should always be greater than the +// numerator. This function checks that. +inline Cost DivideIfGE0(Cost numerator, Cost denominator) { + DCHECK_GE(numerator, -1e-6); + if (numerator < 1e-6) { + return 0.0; + } + return numerator / denominator; +} + +// Dual information related to a SubModel. +// Stores multipliers, reduced costs, and the lower bound, and provides an +// interface that keeps them alligned. +class DualState { + public: + DualState() = default; + template + DualState(const SubModelT& model) + : lower_bound_(.0), + multipliers_(model.num_elements(), .0), + reduced_costs_(model.subset_costs().begin(), + model.subset_costs().end()) {} + + Cost lower_bound() const { return lower_bound_; } + const ElementCostVector& multipliers() const { return multipliers_; } + const SubsetCostVector& reduced_costs() const { return reduced_costs_; } + + // NOTE: This function contains one of the two O(nnz) subgradient steps + template + void DualUpdate(const SubModelT& model, Op multiplier_operator) { + multipliers_.resize(model.num_elements()); + reduced_costs_.resize(model.num_subsets()); + lower_bound_ = .0; + // Update multipliers + for (ElementIndex i : model.ElementRange()) { + multiplier_operator(i, multipliers_[i]); + lower_bound_ += multipliers_[i]; + DCHECK_GE(multipliers_[i], .0); + } + lower_bound_ += ComputeReducedCosts(model, multipliers_, reduced_costs_); + } + + private: + // Single hot point to optimize once for the different use cases. + template + static Cost ComputeReducedCosts(const ModelT& model, + const ElementCostVector& multipliers, + SubsetCostVector& reduced_costs) { + // Compute new reduced costs (O(nnz)) + Cost negative_sum = .0; + for (SubsetIndex j : model.SubsetRange()) { + reduced_costs[j] = model.subset_costs()[j]; + for (ElementIndex i : model.columns()[j]) { + reduced_costs[j] -= multipliers[i]; + } + if (reduced_costs[j] < .0) { + negative_sum += reduced_costs[j]; + } + } + return negative_sum; + } + + Cost lower_bound_; + ElementCostVector multipliers_; + SubsetCostVector reduced_costs_; +}; + +// Utility aggregate to store and pass around both primal and dual states. +struct PrimalDualState { + Solution solution; + DualState dual_state; +}; + +/////////////////////////////////////////////////////////////////////// +///////////////////////////// SUBGRADIENT ///////////////////////////// +/////////////////////////////////////////////////////////////////////// + +// Utilitiy aggregate used by the SubgradientOptimization procedure to +// communicate pass the needed information to the SubgradientCBs interface. +struct SubgradientContext { + const SubModel& model; + const DualState& current_dual_state; + const DualState& best_dual_state; + const Solution& best_solution; + const ElementCostVector& subgradient; +}; + +// Generic set of callbacks hooks used to specialized the behavior of the +// subgradient optimization +class SubgradientCBs { + public: + virtual bool ExitCondition(const SubgradientContext&) = 0; + virtual void RunHeuristic(const SubgradientContext&, Solution&) = 0; + virtual void ComputeMultipliersDelta(const SubgradientContext&, + ElementCostVector& delta_mults) = 0; + virtual bool UpdateCoreModel(SubModel&, PrimalDualState&) = 0; + virtual ~SubgradientCBs() = default; +}; + +// Subgradient optimization procedure. Optimizes the Lagrangian relaxation of +// the Set-Covering problem until a termination criterion si met. +void SubgradientOptimization(SubModel& core_model, SubgradientCBs& cbs, + PrimalDualState& best_state); + +// Subgradient callbacks implementation focused on improving the current best +// dual bound. +class BoundCBs : public SubgradientCBs { + public: + static constexpr Cost kTol = 1e-6; + + BoundCBs(const SubModel& model); + Cost step_size() const { return step_size_; } + bool ExitCondition(const SubgradientContext& context) override; + void ComputeMultipliersDelta(const SubgradientContext& context, + ElementCostVector& delta_mults) override; + void RunHeuristic(const SubgradientContext& context, + Solution& solution) override {} + bool UpdateCoreModel(SubModel& core_model, + PrimalDualState& best_state) override; + + private: + void MakeMinimalCoverageSubgradient(const SubgradientContext& context, + ElementCostVector& subgradient); + + private: + Cost squared_norm_; + + // This addition implements a simplified stabilization technique inspired by + // works in the literature, such as: + // Frangioni, A., Gendron, B., Gorgone, E. (2015): + // "On the computational efficiency of subgradient methods: A case study in + // combinatorial optimization." + // + // The approach aims to reduce oscillations (zig-zagging) in the subgradient + // by using a "moving average" of the current and previous subgradients. The + // current subgradient is weighted by a factor of alpha, while the previous + // subgradients contribution is weighted by (1 - alpha). The parameter alpha + // is set to 0.5 by default but can be adjusted for tuning. The resulting + // stabilized subgradient is referred to as "direction" to distinguish it from + // the original subgradient. + Cost stabilization_coeff = 0.5; // Arbitrary from c4v4 + ElementCostVector direction_; + ElementCostVector prev_direction_; + + std::vector lagrangian_solution_; + + // Stopping condition + Cost prev_best_lb_; + BaseInt max_iter_countdown_; + BaseInt exit_test_countdown_; + BaseInt exit_test_period_; + BaseInt last_core_update_countdown_; + + // Step size + void UpdateStepSize(Cost lower_bound); + Cost step_size_; + Cost last_min_lb_seen_; + Cost last_max_lb_seen_; + BaseInt step_size_update_countdown_; + BaseInt step_size_update_period_; +}; + +//////////////////////////////////////////////////////////////////////// +/////////////////////// MULTIPLIERS BASED GREEDY /////////////////////// +//////////////////////////////////////////////////////////////////////// + +// Higher level function that return a function obtained by the dual multiplier +// based greedy. +Solution RunMultiplierBasedGreedy( + const SubModel& model, const DualState& dual_state, + Cost cost_cutoff = std::numeric_limits::max()); + +// Lower level greedy function which creates or completes a "solution" (seen as +// set of subsets of the current SubModel) until a cutoff size or cost is +// reached. +// Note: since the cutoff might be reached, the returned solution might not be +// feasible. +Cost CoverGreedly(const SubModel& model, const DualState& dual_state, + Cost cost_cutoff, BaseInt size_cutoff, + std::vector& sol_subsets); + +/////////////////////////////////////////////////////////////////////// +//////////////////////// THREE PHASE ALGORITHM //////////////////////// +/////////////////////////////////////////////////////////////////////// + +// Subgradient callbacks implementation focused on wandering near the optimal +// multipliers and invoke the multipliers based greedy heuristic at each +// iteration. +class HeuristicCBs : public SubgradientCBs { + public: + HeuristicCBs() : step_size_(0.1), countdown_(250) {}; + void set_step_size(Cost step_size) { step_size_ = step_size; } + 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(); + 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 { + return false; + } + + private: + Cost step_size_; + BaseInt countdown_; +}; + +PrimalDualState RunThreePhase(SubModel& model, + const Solution& init_solution = {}); + +/////////////////////////////////////////////////////////////////////// +//////////////////////// FULL TO CORE PRICING ///////////////////////// +/////////////////////////////////////////////////////////////////////// + +// 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. +class FullToCoreModel : public SubModel { + using base = SubModel; + struct UpdateTrigger { + BaseInt countdown; + BaseInt period; + BaseInt max_period; + }; + + public: + FullToCoreModel(const Model* full_model); + Cost FixColumns(const std::vector& columns_to_fix) override; + bool UpdateCore(PrimalDualState& core_state) override; + void ResetPricingPeriod(); + const DualState& best_dual_state() const { return best_dual_state_; } + + private: + decltype(auto) IsFocusCol(FullSubsetIndex j) { + return is_focus_col_[static_cast(j)]; + } + decltype(auto) IsFocusRow(FullElementIndex i) { + 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); + + // Views are not composable (for now), so we can either access the full_model + // with the strongly typed view or with the filtered view. + + // 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_); + } + + // Access the full model with the strongly typed view. + StrongModelView StrongTypedFullModelView() const { + return StrongModelView(full_model_); + } + + bool FullToSubModelInvariantCheck(); + + private: + const Model* full_model_; + + // Note: The `is_focus_col_` vector duplicates information already present in + // `SubModelView::cols_sizes_`. However, it does not overlap with any data + // stored in `CoreModel`. Since `CoreModel` is expected to be the primary use + // case, this vector is explicitly maintained here to ensure compatibility. + SubsetBoolVector is_focus_col_; + + // Note: The `is_focus_row_` vector is functionally redundant with either + // `CoreModel::full2core_row_map_` or `SubModelView::rows_sizes_`. These + // existing structures could be used to create the filtered view of the full + // model. However, doing so would require generalizing the current view + // system to work with generic functors instead of vectors of integral types. + // Since the number of elements is assumed to be not prohibitive, a simpler + // implementation that avoids this memory optimization was preferred. + ElementBoolVector is_focus_row_; + + BaseInt num_subsets_; + BaseInt num_elements_; + + DualState full_dual_state_; + DualState best_dual_state_; + + BaseInt update_countdown_; + BaseInt update_period_; + BaseInt update_max_period_; +}; + +// Coverage counter to decide the number of columns to keep in the core model. +static constexpr BaseInt kMinCov = 5; + +} // namespace operations_research::scp + +#endif /* OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H */ diff --git a/ortools/set_cover/set_cover_submodel.cc b/ortools/set_cover/set_cover_submodel.cc new file mode 100644 index 00000000000..02a13eb2559 --- /dev/null +++ b/ortools/set_cover/set_cover_submodel.cc @@ -0,0 +1,299 @@ +// 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/set_cover/set_cover_submodel.h" + +#include "ortools/base/stl_util.h" + +namespace operations_research::scp { + +template +void PrintSubModelSummary(const SubModelT& model) { + std::cout << "SubModelView sizes: rows " << model.num_focus_elements() << "/" + << model.num_elements() << ", columns " << model.num_focus_subsets() + << "/" << model.num_subsets() << '\n'; + std::cout << "Fixing: columns " << model.fixed_columns().size() << " cost " + << model.fixed_cost() << '\n'; +} + +SubModelView::SubModelView(const Model* model) + : base_view(model, &cols_sizes_, &rows_sizes_, &cols_focus_, &rows_focus_), + full_model_(model), + cols_sizes_(model->num_subsets()), + rows_sizes_(model->num_elements()) { + DCHECK(full_model_ != nullptr); + cols_focus_.reserve(model->num_subsets()); + rows_focus_.reserve(model->num_elements()); + for (SubsetIndex j : model->SubsetRange()) { + cols_sizes_[j] = model->columns()[j].size(); + cols_focus_.push_back(j); + } + for (ElementIndex i : model->ElementRange()) { + rows_sizes_[i] = model->rows()[i].size(); + rows_focus_.push_back(i); + } + fixed_columns_.clear(); + fixed_cost_ = 0.0; + PrintSubModelSummary(*this); + DCHECK(ValidateSubModel(*this)); +} + +SubModelView::SubModelView(const Model* model, + const std::vector& columns_focus) + : base_view(model, &cols_sizes_, &rows_sizes_, &cols_focus_, &rows_focus_), + full_model_(model) { + rows_sizes_.resize(full_model_->num_elements(), 0); + for (ElementIndex i : full_model_->ElementRange()) { + rows_sizes_[i] = full_model_->rows()[i].size(); + } + SetFocus(columns_focus); +} + +Cost SubModelView::FixColumns(const std::vector& columns_to_fix) { + DCHECK(full_model_ != nullptr); + Cost old_fixed_cost = fixed_cost_; + if (columns_to_fix.empty()) { + return fixed_cost_ - old_fixed_cost; + } + + for (SubsetIndex j : columns_to_fix) { + DCHECK(cols_sizes_[j] > 0); + fixed_cost_ += full_model_->subset_costs()[j]; + fixed_columns_.push_back(static_cast(j)); + cols_sizes_[j] = 0; + for (ElementIndex i : full_model_->columns()[j]) { + rows_sizes_[i] = 0; + } + } + + gtl::STLEraseAllFromSequenceIf(&cols_focus_, [&](SubsetIndex j) { + DCHECK(full_model_ != nullptr); + if (cols_sizes_[j] > 0) { + cols_sizes_[j] = absl::c_count_if(full_model_->columns()[j], [&](auto i) { + return rows_sizes_[i] > 0; + }); + } + return cols_sizes_[j] == 0; + }); + gtl::STLEraseAllFromSequenceIf( + &rows_focus_, [&](ElementIndex i) { return !rows_sizes_[i]; }); + + PrintSubModelSummary(*this); + DCHECK(ValidateSubModel(*this)); + return fixed_cost_ - old_fixed_cost; +} + +void SubModelView::SetFocus(const std::vector& columns_focus) { + DCHECK(full_model_ != nullptr); + DCHECK(!rows_sizes_.empty()); + if (columns_focus.empty()) { + return; + } + cols_focus_.clear(); + rows_focus_.clear(); + + ElementBoolVector enabled_rows(full_model_->num_elements(), false); + for (ElementIndex i : full_model_->ElementRange()) { + enabled_rows[i] = rows_sizes_[i] > 0; + } + cols_sizes_.assign(full_model_->num_subsets(), 0); + rows_sizes_.assign(full_model_->num_elements(), 0); + for (FullSubsetIndex full_j : columns_focus) { + SubsetIndex j = static_cast(full_j); + for (ElementIndex i : full_model_->columns()[j]) { + if (enabled_rows[i] > 0) { + ++cols_sizes_[j]; + ++rows_sizes_[i]; + } + } + if (cols_sizes_[j] > 0) { + cols_focus_.push_back(j); + } + } + for (ElementIndex i : full_model_->ElementRange()) { + if (rows_sizes_[i] > 0) { + 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()) { + 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()); +} + +CoreModel::CoreModel(const Model* model, + const std::vector& columns_focus) + : Model(), + full_model_(model), + core2full_row_map_(model->num_elements()), + full2core_row_map_(model->num_elements()) { + 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()); + SetFocus(columns_focus); +} + +// 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 +// using "cleanish" code. +void CoreModel::SetFocus(const std::vector& columns_focus) { + if (columns_focus.empty()) { + return; + } + + // TODO(c4v4): change model in-place to avoid reallocations. + Model& submodel = static_cast(*this); + submodel = Model(); + core2full_col_map_.clear(); + + // 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]; + if (core_i != null_element_index) { + if (first_row) { + // SetCoverModel lacks a way to remove columns + first_row = false; + submodel.AddEmptySubset(full_model_.subset_costs()[full_j]); + } + submodel.AddElementToLastSubset(core_i); + } + } + } + + submodel.CreateSparseRowView(); + PrintSubModelSummary(*this); + DCHECK(ValidateSubModel(*this)); +} + +// Mark columns and row that will be removed from the core model +// The "to-be-removed" indices are marked by setting the relative core->full +// mappings to null_*_index. +void CoreModel::MarkNewFixingInMaps( + const std::vector& columns_to_fix) { + for (SubsetIndex old_core_j : columns_to_fix) { + fixed_cost_ += subset_costs()[old_core_j]; + fixed_columns_.push_back(core2full_col_map_[old_core_j]); + + core2full_col_map_[old_core_j] = null_full_subset_index; + for (ElementIndex old_core_i : columns()[old_core_j]) { + core2full_row_map_[old_core_i] = null_full_element_index; + } + } +} + +// Once fixed columns and covered rows are marked, we need to create a new row +// mapping, both core->full(returned) and full->core(modifed in-place) +CoreToFullElementMapVector CoreModel::MakeOrFillBothRowMaps() { + full2core_row_map_.assign(full_model_.num_elements(), null_element_index); + CoreToFullElementMapVector new_c2f_row_map; // Future core->full mapping. + for (ElementIndex old_core_i : ElementRange()) { + FullElementIndex full_i = core2full_row_map_[old_core_i]; + if (full_i != null_full_element_index) { + full2core_row_map_[full_i] = ElementIndex(new_c2f_row_map.size()); + new_c2f_row_map.push_back(full_i); + } + } + return new_c2f_row_map; +} + +// Create a new core model by applying the remapping from the old core model to +// the new one considering the given column fixing. +// Both the old and new core->full row mappings are required too keep track of +// what changed, the old mapping gets overwritten with the new one at the end. +// Empty columns are detected and removed - or better - not added. +Model CoreModel::MakeNewCoreModel( + const CoreToFullElementMapVector& new_c2f_row_map) { + Model new_submodel; + BaseInt new_core_j = 0; + // Loop over old core column indices. + for (SubsetIndex old_core_j : SubsetRange()) { + // If the column is not marked, then it should be mapped. + FullSubsetIndex full_j = core2full_col_map_[old_core_j]; + if (full_j != null_full_subset_index) { + bool first_row = true; + // Loop over the old core column (with old core row indices). + for (ElementIndex old_core_i : columns()[old_core_j]) { + // If the row is not marked, then it should be mapped. + FullElementIndex full_i = core2full_row_map_[old_core_i]; + if (full_i != null_full_element_index) { + if (first_row) { + // SetCoverModel lacks a way to remove columns + first_row = false; + new_submodel.AddEmptySubset(full_model_.subset_costs()[full_j]); + + // Put the full index in the proper (new) position. + // Note that old_core_j >= new_core_j is always true. + SubsetIndex new_j(new_core_j++); + core2full_col_map_[new_j] = full_j; + } + ElementIndex new_core_i = full2core_row_map_[full_i]; + DCHECK(new_core_i != null_element_index); + new_submodel.AddElementToLastSubset(new_core_i); + } + } + } + } + + core2full_col_map_.resize(new_core_j); + core2full_row_map_ = std::move(new_c2f_row_map); + new_submodel.CreateSparseRowView(); + + return new_submodel; +} + +Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { + if (columns_to_fix.empty()) { + return .0; + } + Cost old_fixed_cost = fixed_cost_; + + // Mark columns to be fixed and rows that will be covered by them + MarkNewFixingInMaps(columns_to_fix); + + // Compute new core->full(returned) and full->core(modified original) row maps + CoreToFullElementMapVector new_c2f_row_map = MakeOrFillBothRowMaps(); + + // 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_)); + + return fixed_cost_ - old_fixed_cost; +} + +} // 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 new file mode 100644 index 00000000000..18dc20d29e9 --- /dev/null +++ b/ortools/set_cover/set_cover_submodel.h @@ -0,0 +1,289 @@ +// 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_SET_COVER_SET_COVER_SUBMODEL_H +#define OR_TOOLS_SET_COVER_SET_COVER_SUBMODEL_H + +#include "ortools/set_cover/set_cover_views.h" + +namespace operations_research::scp { + +// TODO(anyone): since we are working within the scp namespace, the "SetCover*" +// prefix became redundant and can be removed. For now, only redefine it +// locally. +using Model = SetCoverModel; + +// Forward declarations, see below for the definition of the classes. +struct PrimalDualState; +struct Solution; + +// The CFT algorithm generates sub-models in two distinct ways: +// +// 1. It fixes specific columns (incrementally) into any generated solution. +// Once a column is fixed, it is excluded from future decisions, as it is +// already part of the solution. Additionally, rows that are covered by the +// fixed columns are removed from consideration as well, along with any +// columns that exclusively cover those rows, as they become redundant. The +// fixing process starts with the entire model and progressively fixes more +// columns until it becomes empty. A "view-based" sub-model is well-suited +// for this part. +// +// 2. The CFT mostly works on a "core" sub-model by focusing on a subset of +// columns. The core model is derived from the original model but is +// significantly smaller, as it typically includes only a limited number of +// columns per row (on average, around six columns per row). Unlike the +// incremental nature of column fixing, core models are constructed from +// scratch during each update. This type of small model can take advantage of +// a Model object which stores the sub-model explicitly in memory, avoiding +// looping over "inactive" columns and rows. Both SubModelView and CoreModel +// can be used as a core model. +// +// Two types of "core-model" representations are implemented, both of which can +// be used interchangeably: +// +// 1. SubModelView: A lightweight view of the original model. It dynamically +// filters and exposes only the active rows and columns from the original +// data structures, skipping "inactive" items. +// +// 2. CoreModel: A fully compacted and explicit representation of a sub-model. +// It stores the filtered data explicitly, making it more suitable +// for scenarios where compact storage and faster access are required. +// +// While CoreModel stores an explicit representation of the sub-model, +// SubModelView maintains vectors sized according to the original model's +// dimensions. As a result, depending on the dimensions of the original model, +// CoreModel can actually be more memory-efficient. +class SubModelView; +class CoreModel; +using SubModel = CoreModel; + +// `SubModelView` provides a mechanism to interact with a subset of the rows and +// columns of a SetCoverModel, effectively creating a filtered view of the +// model. This abstraction allows operations to be performed on a restricted +// portion of the model without modifying the original data structure. The +// filtering is achieved using index lists and sizes vectors, which define the +// active rows and columns. This approach ensures flexibility and avoids +// unnecessary duplication of data. Columns/rows sizes are uses to both keep +// track of the number of elements in them and also provide the "activation" +// status: (item size == 0) <==> inactive +// SubModelView inherits from IndexListSubModelView, which provides the "view" +// machinery. +class SubModelView : public IndexListModelView { + using base_view = IndexListModelView; + + public: + // Empty initialization to facilitate delayed construction + SubModelView() = default; + + // Identity sub-model: all items are considered + SubModelView(const Model* model); + + // Focus construction: create a sub-model with only the required items + SubModelView(const Model* model, + const std::vector& columns_focus); + + virtual ~SubModelView() = default; + + ///////// Core-model interface: ///////// + + // Current fixed cost: sum of the cost of the fixed columns + Cost fixed_cost() const { return fixed_cost_; } + + // List of fixed columns. + const std::vector& fixed_columns() const { + return fixed_columns_; + } + + // Redefine the active items. The new sub-model will ignore all columns not in + // focus and (optionally) the rows for which row_flags is not true. It does + // not overwrite the current fixing. + void SetFocus(const std::vector& columns_focus); + + // 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); + + // 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; } + + private: + // Pointer to the original model + const Model* full_model_; + + // Columns/rows sizes after filtering (size==0 <==> inactive) + SubsetToIntVector cols_sizes_; + ElementToIntVector rows_sizes_; + + // List of columns/rows currectly active + std::vector cols_focus_; + std::vector rows_focus_; + + // Fixing data + std::vector fixed_columns_; + Cost fixed_cost_ = .0; +}; + +// CoreModel stores a subset of the filtered columns and rows in an explicit +// Model object. +// The indices are compacted and mapped to the range [0, ], +// effectively creating a smaller set-covering model. Similar to SubModelView, +// the core model supports column fixing and focusing on a subset of the +// original model. Mappings are maintained to translate indices back to the +// original model space. +class CoreModel : private Model { + public: + // Empty initialization to facilitate delayed construction + CoreModel() = default; + + // Identity sub-model: all items are considered + CoreModel(const Model* model); + + // Focus construction: create a sub-model with only the required items + CoreModel(const Model* model, + const std::vector& columns_focus); + + virtual ~CoreModel() = default; + + ///////// Sub-model view interface: ///////// + BaseInt num_subsets() const { return full_model_.num_subsets(); } + BaseInt num_elements() const { return full_model_.num_elements(); } + BaseInt num_focus_subsets() const { return Model::num_subsets(); } + BaseInt num_focus_elements() const { return Model::num_elements(); } + BaseInt column_size(SubsetIndex j) const { + DCHECK(SubsetIndex() <= j && j < SubsetIndex(num_subsets())); + return columns()[j].size(); + } + BaseInt row_size(ElementIndex i) const { + DCHECK(ElementIndex() <= i && i < ElementIndex(num_elements())); + return rows()[i].size(); + } + + FullElementIndex MapCoreToFullElementIndex(ElementIndex core_i) const { + DCHECK(ElementIndex() <= core_i && core_i < ElementIndex(num_elements())); + return core2full_row_map_[core_i]; + } + 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 { + DCHECK(SubsetIndex() <= core_j && core_j < SubsetIndex(num_subsets())); + return core2full_col_map_[core_j]; + } + // Member function relevant for the CFT inherited from Model + using Model::columns; + using Model::ElementRange; + using Model::rows; + using Model::subset_costs; + using Model::SubsetRange; + + ///////// Core-model interface: ///////// + + // Current fixed cost: sum of the cost of the fixed columns + Cost fixed_cost() const { return fixed_cost_; } + // List of fixed columns. + const std::vector& fixed_columns() const { + return fixed_columns_; + } + + // Redefine the active items. The new sub-model will ignore all columns not in + // focus and (optionally) the rows for which row_flags is not true. It does + // not overwrite the current fixing. + void SetFocus(const std::vector& columns_focus); + + // 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); + + // 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; } + + private: + void MarkNewFixingInMaps(const std::vector& columns_to_fix); + CoreToFullElementMapVector MakeOrFillBothRowMaps(); + Model MakeNewCoreModel(const CoreToFullElementMapVector& new_c2f_col_map); + + // Pointer to the original model + StrongModelView full_model_; + + FullToCoreElementMapVector full2core_row_map_; + CoreToFullElementMapVector core2full_row_map_; + CoreToFullSubsetMapVector core2full_col_map_; + + // Fixing data + Cost fixed_cost_ = .0; + std::vector fixed_columns_; + + static constexpr SubsetIndex null_subset_index = + std::numeric_limits::max(); + static constexpr ElementIndex null_element_index = + std::numeric_limits::max(); + static constexpr FullSubsetIndex null_full_subset_index = + std::numeric_limits::max(); + static constexpr FullElementIndex null_full_element_index = + std::numeric_limits::max(); +}; + +template +bool ValidateSubModel(const SubModelT& model) { + if (model.num_elements() <= 0) { + std::cerr << "Sub-Model has no elements.\n"; + return false; + } + if (model.num_subsets() <= 0) { + std::cerr << "Sub-Model has no subsets.\n"; + return false; + } + + for (SubsetIndex j : model.SubsetRange()) { + const auto& column = model.columns()[j]; + if (model.column_size(j) == 0) { + std::cerr << "Column " << j << " is empty.\n"; + return false; + } + BaseInt j_size = std::distance(column.begin(), column.end()); + if (j_size != model.column_size(j)) { + std::cerr << "Sub-Model size mismatch on column " << j << ", " << j_size + << " != " << model.column_size(j) << "\n"; + return false; + } + } + + for (ElementIndex i : model.ElementRange()) { + const auto& row = model.rows()[i]; + if (model.row_size(i) == 0) { + std::cerr << "Row " << i << " is empty.\n"; + return false; + } + BaseInt i_size = std::distance(row.begin(), row.end()); + if (i_size != model.row_size(i)) { + std::cerr << "Sub-Model size mismatch on row " << i << ", " << i_size + << " != " << model.row_size(i) << "\n"; + return false; + } + } + return true; +} + +} // namespace operations_research::scp +#endif /* OR_TOOLS_SET_COVER_SET_COVER_SUBMODEL_H */ diff --git a/ortools/set_cover/set_cover_views.h b/ortools/set_cover/set_cover_views.h new file mode 100644 index 00000000000..d9c97d725c1 --- /dev/null +++ b/ortools/set_cover/set_cover_views.h @@ -0,0 +1,247 @@ +// 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_SET_COVER_SET_COVER_VIEWS_H +#define OR_TOOLS_SET_COVER_SET_COVER_VIEWS_H + +#include + +#include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_model.h" +#include "views.h" + +namespace operations_research { + +// In the CFT algorithm, indices from different models are frequently used, +// and mixing them can lead to errors. To prevent such mistakes, strong-typed +// wrappers are employed. There are three available approaches for handling +// these indices: +// 1. Full-model strong-typed indices + {Subset,Element}Index for the core model +// 2. Core-model strong-typed indices + {Subset,Element}Index for the full model +// 3. Define new strong-typed indices both full-model and core-model +// Introducing a new set of strong-typed indices, however, can lead to a +// cascade of code duplication (or template proliferation). It also requires +// additional "view" boilerplate to properly handle the different types, +// increasing complexity. +// Currently, the simplest approach is to define only full-model indices while +// reusing the original strong types for the core model. The main challenge +// arises in FullToCoreModel, where a "filtered" full-model must be handled. In +// such cases, static casts are employed to manage the type conversions +// effectively. +DEFINE_STRONG_INT_TYPE(FullSubsetIndex, BaseInt); +DEFINE_STRONG_INT_TYPE(FullElementIndex, BaseInt); + +// Syntactic sugar to define strong typed indices casts. +// Note: look at `strong_int.h` for more details about `StrongIntConvert` +#define ENABLE_EXPLICIT_STRONG_TYPE_CAST(FROM, TO) \ + constexpr TO StrongIntConvert(FROM j, TO* /*unused*/) { \ + return TO(static_cast(j)); \ + } + +ENABLE_EXPLICIT_STRONG_TYPE_CAST(SubsetIndex, FullSubsetIndex); +ENABLE_EXPLICIT_STRONG_TYPE_CAST(FullSubsetIndex, SubsetIndex); +ENABLE_EXPLICIT_STRONG_TYPE_CAST(ElementIndex, FullElementIndex); +ENABLE_EXPLICIT_STRONG_TYPE_CAST(FullElementIndex, ElementIndex); + +using FullElementCostVector = util_intops::StrongVector; +using FullSubsetCostVector = util_intops::StrongVector; +using FullElementBoolVector = util_intops::StrongVector; +using FullSubsetBoolVector = 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 +// ElementIndex mappings are stored to translate back to the original model +// space. +using FullToCoreElementMapVector = + util_intops::StrongVector; +using CoreToFullElementMapVector = + util_intops::StrongVector; + +// The same applies to SubsetIndex, which also needs to be mapped back to the +// original indexing space. +using FullToCoreSubsetMapVector = + util_intops::StrongVector; +using CoreToFullSubsetMapVector = + util_intops::StrongVector; + +class StrongModelView { + private: + // Transformations to convert between the core and full model columns. + struct SparseColTransform { + auto operator()(const SparseColumn& column) const + -> util_intops::TransformView< + ElementIndex, ColumnEntryIndex, + util_intops::TypeCastTransform> { + return {&column}; + } + }; + + // Transformations to convert between the core and full model rows. + struct SparseRowTransform { + auto operator()(const SparseRow& row) const -> util_intops::TransformView< + SubsetIndex, RowEntryIndex, + util_intops::TypeCastTransform> { + return {&row}; + } + }; + + public: + StrongModelView() = default; + StrongModelView(const SetCoverModel* model) : model_(model) {} + + BaseInt num_subsets() const { return model_->num_subsets(); } + BaseInt num_elements() const { return model_->num_elements(); } + + auto subset_costs() const + -> util_intops::TransformView { + return {&model_->subset_costs()}; + } + auto columns() const + -> util_intops::TransformView { + return {&model_->columns()}; + } + auto rows() const -> util_intops::TransformView { + return {&model_->rows()}; + } + auto SubsetRange() const -> util_intops::StrongIntRange { + return {FullSubsetIndex(), FullSubsetIndex(num_subsets())}; + } + auto ElementRange() const -> util_intops::StrongIntRange { + return {FullElementIndex(), FullElementIndex(num_elements())}; + } + + private: + const SetCoverModel* model_; +}; + +class IndexListModelView { + public: + IndexListModelView() = default; + IndexListModelView(const SetCoverModel* model, + const SubsetToIntVector* cols_sizes, + const ElementToIntVector* rows_sizes, + const std::vector* cols_focus, + const std::vector* rows_focus) + : model_(model), + cols_sizes_(cols_sizes), + rows_sizes_(rows_sizes), + cols_focus_(cols_focus), + rows_focus_(rows_focus) {} + + BaseInt num_subsets() const { return model_->num_subsets(); } + BaseInt num_elements() const { return model_->num_elements(); } + BaseInt num_focus_subsets() const { return cols_focus_->size(); } + BaseInt num_focus_elements() const { return rows_focus_->size(); } + + auto subset_costs() const -> util_intops::IndexListView { + return {&model_->subset_costs(), cols_focus_}; + } + auto columns() const -> util_intops::TwoLevelsView< + util_intops::IndexListView, + ElementToIntVector> { + return {{&model_->columns(), cols_focus_}, rows_sizes_}; + } + auto rows() const -> util_intops::TwoLevelsView< + util_intops::IndexListView, SubsetToIntVector> { + return {{&model_->rows(), rows_focus_}, cols_sizes_}; + } + const std::vector& SubsetRange() const { return *cols_focus_; } + const std::vector& ElementRange() const { return *rows_focus_; } + FullElementIndex MapCoreToFullElementIndex(ElementIndex core_i) const { + DCHECK(ElementIndex() <= core_i && core_i < ElementIndex(num_elements())); + return static_cast(core_i); + } + ElementIndex MapFullToCoreElementIndex(FullElementIndex full_i) const { + DCHECK(FullElementIndex() <= full_i && + full_i < FullElementIndex(num_elements())); + return static_cast(full_i); + } + FullSubsetIndex MapCoreToFullSubsetIndex(SubsetIndex core_j) const { + DCHECK(SubsetIndex() <= core_j && core_j < SubsetIndex(num_subsets())); + return static_cast(core_j); + } + BaseInt column_size(SubsetIndex j) const { + DCHECK(SubsetIndex() <= j && j < SubsetIndex(num_subsets())); + return (*cols_sizes_)[j]; + } + BaseInt row_size(ElementIndex i) const { + DCHECK(ElementIndex() <= i && i < ElementIndex(num_elements())); + return (*rows_sizes_)[i]; + } + + private: + const SetCoverModel* model_; + const SubsetToIntVector* cols_sizes_; + const ElementToIntVector* rows_sizes_; + const std::vector* cols_focus_; + const std::vector* rows_focus_; +}; + +// A lightweight sub-model view that uses boolean vectors to enable or disable +// specific items. Iterating over all active columns or rows is less efficient, +// particularly when only a small subset is active. +// NOTE: this view does **not** store any size-related information. +class FilterModelView { + public: + FilterModelView() = default; + FilterModelView(const SetCoverModel* model, + const SubsetBoolVector* cols_sizes, + const ElementBoolVector* rows_sizes, BaseInt num_subsets, + BaseInt num_elements) + : model_(model), + is_focus_col_(cols_sizes), + is_focus_row_(rows_sizes), + num_subsets_(num_subsets), + num_elements_(num_elements) {} + + BaseInt num_subsets() const { return model_->num_subsets(); } + BaseInt num_elements() const { return model_->num_elements(); } + BaseInt num_focus_subsets() const { return num_subsets_; } + BaseInt num_focus_elements() const { return num_elements_; } + + auto subset_costs() const + -> util_intops::IndexFilterView { + return {&model_->subset_costs(), is_focus_col_}; + } + auto columns() const -> util_intops::TwoLevelsView< + util_intops::IndexFilterView, + ElementBoolVector> { + return {{&model_->columns(), is_focus_col_}, is_focus_row_}; + } + auto rows() const -> util_intops::TwoLevelsView< + util_intops::IndexFilterView, + SubsetBoolVector> { + return {{&model_->rows(), is_focus_row_}, is_focus_col_}; + } + auto SubsetRange() const + -> util_intops::FilterIndexRangeView { + return {is_focus_col_}; + } + auto ElementRange() const + -> util_intops::FilterIndexRangeView { + return {is_focus_row_}; + } + + private: + const SetCoverModel* model_; + const SubsetBoolVector* is_focus_col_; + const ElementBoolVector* is_focus_row_; + BaseInt num_subsets_; + BaseInt num_elements_; +}; + +} // namespace operations_research + +#endif /* OR_TOOLS_SET_COVER_SET_COVER_VIEWS_H */ diff --git a/ortools/set_cover/views.h b/ortools/set_cover/views.h new file mode 100644 index 00000000000..ba4c5a46c38 --- /dev/null +++ b/ortools/set_cover/views.h @@ -0,0 +1,470 @@ +// 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_SET_COVER_VIEWS_H +#define OR_TOOLS_ORTOOLS_SET_COVER_VIEWS_H + +#include +#include + +#include + +// 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 +// Cover problem, where we often work with subsets of columns or rows. Despite +// this, each column or row still contains the original indices, which are +// used for referring to other columns/rows. +// +// This design abstracts access to the underlying container, ensuring consistent +// behavior across the following scenarios: +// 1. Directly using the original container. +// 2. Accessing a subset of the original items via a view. +// 3. Using a new container with a compacted subset of items, indexing them with +// the position the have in the new container and not in the original one. +// This also need the old indices stored in rows/columns to be mapped into +// the new ones. + +namespace util_intops { + +namespace util { + +template +using range_const_iterator_type = + absl::remove_cvref_t().begin())>; +template +using range_value_type = absl::remove_cvref_t< + decltype(*std::declval>())>; + +// Enable compatibility with STL algorithms. +template +struct IteratorCRTP { + using iterator_category = std::input_iterator_tag; + using value_type = ValueT; + using difference_type = std::ptrdiff_t; + using pointer = IterT; + using reference = value_type&; + bool operator==(const IterT& other) const { + return !(*static_cast(this) != other); + } + pointer operator->() const { return static_cast(this); } +}; + +template +decltype(auto) at(const ValueRangeT* container, IndexT index) { + DCHECK(IndexT(0) <= index && index < IndexT(container->size())); + return (*container)[index]; +} + +} // namespace util + +// View exposing only the indices of a container that are marked as active. +// Looping over this view is equivalent to: +// +// for (IntT index; index < IntTsize()); ++index) { +// if (is_active[index]) { +// your_code(index); +// } +// } +// +template +class FilterIndexRangeView { + public: + struct FilterIndicesViewIterator + : util::IteratorCRTP { + FilterIndicesViewIterator(IntT index, const EnableVectorT* is_active) + : index_(index), is_active_(is_active) { + AdjustToValidValue(); + std::vector vb; + } + bool operator!=(FilterIndicesViewIterator other) const { + return index_ != other.index_; + } + FilterIndicesViewIterator& operator++() { + ++index_; + AdjustToValidValue(); + return *this; + } + IntT operator*() const { return index_; } + + private: + void AdjustToValidValue() { + while (index_ < IntT(is_active_->size()) && + !util::at(is_active_, index_)) { + ++index_; + } + } + + IntT index_; + const EnableVectorT* is_active_; + }; + + FilterIndexRangeView(const EnableVectorT* is_active) + : is_active_(is_active) {} + FilterIndicesViewIterator begin() const { + return FilterIndicesViewIterator(IntT(0), is_active_); + } + FilterIndicesViewIterator end() const { + return FilterIndicesViewIterator(IntT(is_active_->size()), is_active_); + } + + private: + const EnableVectorT* is_active_; +}; + +// View exposing only the elements of a container that are indexed by a list of +// indices. +// Looping over this view is equivalent to: +// +// for (decltype(auto) index : indices) { +// your_code(container[index]); +// } +// +template +class IndexListView { + public: + using value_type = const ValueT; + using index_type = const IndexT; + using value_iterator = typename absl::Span::iterator; + using index_iterator = typename absl::Span::iterator; + + struct IndexListViewIterator + : util::IteratorCRTP { + IndexListViewIterator(absl::Span values, index_iterator iter) + : values_(values), index_iter_(iter) {} + bool operator!=(const IndexListViewIterator& other) const { + DCHECK_EQ(values_.data(), other.values_.data()); + return index_iter_ != other.index_iter_; + } + IndexListViewIterator& operator++() { + ++index_iter_; + return *this; + } + decltype(auto) operator*() const { + return values_[static_cast(*index_iter_)]; + } + + private: + absl::Span values_; + index_iterator index_iter_; + }; + + IndexListView() = default; + + template + IndexListView(const ValueRangeT* values, const IndexRangeT* indices) + : values_(absl::MakeConstSpan(*values)), + indices_(absl::MakeConstSpan(*indices)) {} + + auto size() const { return indices_.size(); } + bool empty() const { return indices_.empty(); } + + // NOTE: uses indices of the original container, not the filtered one + decltype(auto) operator[](index_type index) const { + // TODO(user): we could check that index is in indices_, but that's is O(n). + return values_[static_cast(index)]; + } + IndexListViewIterator begin() const { + return IndexListViewIterator(values_, indices_.begin()); + } + IndexListViewIterator end() const { + return IndexListViewIterator(values_, indices_.end()); + } + + private: + absl::Span values_; + absl::Span indices_; +}; + +// This view provides access to elements of a container of integral types, which +// are used as indices to a filter (in) vector. +// The view filters and returns only those elements that, when used a subscript +// to the filter vector, are evaluated to a true value. +// Looping over this view is equivalent to: +// +// for (decltype(auto) item : container) { +// if (is_active[item]) { +// your_code(iter); +// } +// } +// +template +class ValueFilterView { + public: + using value_type = const ValueT; + using value_iterator = typename absl::Span::iterator; + + struct ValueFilterViewIterator + : util::IteratorCRTP { + ValueFilterViewIterator(value_iterator iterator, value_iterator end, + const EnableVectorT* is_active) + : iterator_(iterator), end_(end), is_active_(is_active) { + AdjustToValidValue(); + } + bool operator!=(const ValueFilterViewIterator& other) const { + DCHECK_EQ(is_active_, other.is_active_); + return iterator_ != other.iterator_; + } + ValueFilterViewIterator& operator++() { + ++iterator_; + AdjustToValidValue(); + return *this; + } + decltype(auto) operator*() const { return *iterator_; } + + private: + void AdjustToValidValue() { + while (iterator_ != end_ && !util::at(is_active_, *iterator_)) { + ++iterator_; + } + } + + value_iterator iterator_; + value_iterator end_; + const EnableVectorT* is_active_; + }; + + template + ValueFilterView(const ValueRangeT* values, const EnableVectorT* is_active) + : values_(absl::MakeConstSpan(*values)), is_active_(is_active) { + DCHECK(values != nullptr); + DCHECK(is_active != nullptr); + } + ValueFilterViewIterator begin() const { + return ValueFilterViewIterator(values_.begin(), values_.end(), is_active_); + } + ValueFilterViewIterator end() const { + return ValueFilterViewIterator(values_.end(), values_.end(), is_active_); + } + + // NOTE: uses indices of the original container, not the filtered one + template + decltype(auto) operator[](IndexT index) const { + decltype(auto) value = values_[static_cast(index)]; + DCHECK(util::at(is_active_, value)) + << "Inactive value. Are you using relative indices?"; + return value; + } + + private: + absl::Span values_; + const EnableVectorT* is_active_; +}; + +// Somewhat equivalent to ValueFilterView +// Looping over this view is equivalent to: +// +// auto c_it = container.begin(); +// auto active_it = is_active.begin(); +// for (; active_it != is_active.end(); ++active_it, ++c_it) { +// if (*active_it) { +// your_code(*c_it); +// } +// } +// +template +class IndexFilterView { + public: + using value_type = const ValueT; + using value_iterator = typename absl::Span::iterator; + using enable_iterator = util::range_const_iterator_type; + + struct IndexFilterViewIterator + : util::IteratorCRTP { + IndexFilterViewIterator(value_iterator iterator, + enable_iterator is_active_begin, + enable_iterator is_active_end) + : iterator_(iterator), + is_active_iter_(is_active_begin), + is_active_end_(is_active_end) { + AdjustToValidValue(); + } + bool operator!=(const IndexFilterViewIterator& other) const { + DCHECK(is_active_end_ == other.is_active_end_); + return iterator_ != other.iterator_; + } + IndexFilterViewIterator& operator++() { + ++is_active_iter_; + ++iterator_; + AdjustToValidValue(); + return *this; + } + decltype(auto) operator*() const { return *iterator_; } + + private: + void AdjustToValidValue() { + while (is_active_iter_ != is_active_end_ && !*is_active_iter_) { + ++is_active_iter_; + ++iterator_; + } + } + + value_iterator iterator_; + enable_iterator is_active_iter_; + enable_iterator is_active_end_; + }; + + template + IndexFilterView(const ValueRangeT* values, const EnableVectorT* is_active_) + : values_(absl::MakeConstSpan(*values)), is_active_(is_active_) { + DCHECK(values != nullptr); + DCHECK(is_active_ != nullptr); + DCHECK_EQ(values->size(), is_active_->size()); + } + IndexFilterViewIterator begin() const { + return IndexFilterViewIterator(values_.begin(), is_active_->begin(), + is_active_->end()); + } + IndexFilterViewIterator end() const { + return IndexFilterViewIterator(values_.end(), is_active_->end(), + is_active_->end()); + } + + // NOTE: uses indices of the original container, not the filtered one + template + decltype(auto) operator[](IndexT index) const { + DCHECK(util::at(is_active_, index)) + << "Inactive value. Are you using relative indices?"; + return values_[static_cast(index)]; + } + + private: + absl::Span values_; + const EnableVectorT* is_active_; +}; + +// This view provides a mechanism to access and filter elements in a 2D +// container. The filtering is applied in two stages: +// 1. The first dimension is generic and can be both and index-list or +// bool-vector based view. +// 2. The second dimension (items of each sub-container) is further filtered +// using a boolean-vector-like object, which determines which elements within +// the sub-container are included in the view. +template +class TwoLevelsView : public Lvl1ViewT { + public: + using level1_iterator = util::range_const_iterator_type; + using level1_value = util::range_value_type; + using level2_value = util::range_value_type; + using level2_type = ValueFilterView; + + struct TwoLevelsViewIterator + : util::IteratorCRTP { + TwoLevelsViewIterator(level1_iterator iter, + const EnableVectorT* active_items) + : iter_(iter), active_items_(active_items) {} + bool operator!=(const TwoLevelsViewIterator& other) const { + DCHECK_EQ(active_items_, other.active_items_); + return iter_ != other.iter_; + } + TwoLevelsViewIterator& operator++() { + ++iter_; + return *this; + } + level2_type operator*() const { + return level2_type(&(*iter_), active_items_); + } + + private: + level1_iterator iter_; + const EnableVectorT* active_items_; + }; + + TwoLevelsView() = default; + TwoLevelsView(Lvl1ViewT lvl1_view, const EnableVectorT* active_items) + : Lvl1ViewT(lvl1_view), active_items_(active_items) {} + TwoLevelsViewIterator begin() const { + return TwoLevelsViewIterator(Lvl1ViewT::begin(), active_items_); + } + TwoLevelsViewIterator end() const { + return TwoLevelsViewIterator(Lvl1ViewT::end(), active_items_); + } + + template + level2_type operator[](indexT i) const { + auto& level2_container = Lvl1ViewT::operator[](i); + return level2_type(&level2_container, active_items_); + } + + private: + const EnableVectorT* active_items_; +}; + +struct NoTransform { + template + T&& operator()(T&& value) const { + return std::forward(value); + } +}; + +template +struct TypeCastTransform { + ToT operator()(FromT value) const { return static_cast(value); } +}; + +// View applying stateless transformation to the values of a continugous +// container. Looping over this view is equivalent to: +// +// for (IndexT i(0); i < IndexT(container.size()); ++i) { +// your_code(ValueTransformT()(values_[static_cast(i)])); +// } +// +template +class TransformView { + public: + using value_type = const ValueT; + using value_iterator = typename absl::Span::iterator; + + struct TransformViewIterator + : util::IteratorCRTP { + TransformViewIterator(value_iterator iterator) : iterator_(iterator) {} + bool operator!=(const TransformViewIterator& other) const { + return iterator_ != other.iterator_; + } + TransformViewIterator& operator++() { + ++iterator_; + return *this; + } + decltype(auto) operator*() const { return ValueTransformT()(*iterator_); } + + private: + value_iterator iterator_; + }; + + TransformView() = default; + + template + TransformView(const ValueRangeT* values) + : values_(absl::MakeConstSpan(*values)) {} + + auto size() const { return values_.size(); } + bool empty() const { return values_.empty(); } + + decltype(auto) operator[](IndexT index) const { + return ValueTransformT()(values_[static_cast(index)]); + } + TransformViewIterator begin() const { + return TransformViewIterator(values_.begin()); + } + TransformViewIterator end() const { + return TransformViewIterator(values_.end()); + } + + private: + absl::Span values_; +}; + +} // namespace util_intops + +#endif /* OR_TOOLS_ORTOOLS_SET_COVER_VIEWS_H */