diff --git a/ortools/algorithms/bin_packing.cc b/ortools/algorithms/bin_packing.cc new file mode 100644 index 00000000000..0f512b131b9 --- /dev/null +++ b/ortools/algorithms/bin_packing.cc @@ -0,0 +1,501 @@ + +// Copyright 2025 Francesco Cavaliere +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/algorithms/bin_packing.h" + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "ortools/base/stl_util.h" +#include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_cft.h" +#include "ortools/set_cover/set_cover_submodel.h" +#include "ortools/util/filelineiter.h" + +namespace operations_research { + +void BinPackingModel::set_bin_capacity(Cost capacity) { + if (capacity <= 0) { + LOG(WARNING) << "Bin capacity must be positive."; + return; + } + bin_capcaity_ = capacity; +} +void BinPackingModel::AddItem(Cost weight) { + if (weight > bin_capcaity_) { + LOG(WARNING) << "Element weight exceeds bin capacity."; + return; + } + weigths_.push_back(weight); + is_sorted_ = false; +} +void BinPackingModel::SortWeights() { + if (!is_sorted_) { + absl::c_sort(weigths_); + is_sorted_ = true; + } +} + +namespace { +template +bool SimpleAtoValue(absl::string_view str, int_type* out) { + return absl::SimpleAtoi(str, out); +} +bool SimpleAtoValue(absl::string_view str, bool* out) { + return absl::SimpleAtob(str, out); +} +bool SimpleAtoValue(absl::string_view str, double* out) { + return absl::SimpleAtod(str, out); +} +bool SimpleAtoValue(absl::string_view str, float* out) { + return absl::SimpleAtof(str, out); +} +} // namespace + +BinPackingModel ReadBpp(absl::string_view filename) { + BinPackingModel model; + BaseInt num_items = 0; + for (const std::string& line : + FileLines(filename, FileLineIterator::REMOVE_INLINE_CR | + FileLineIterator::REMOVE_BLANK_LINES)) { + if (num_items == 0) { + if (!SimpleAtoValue(line, &num_items)) { + LOG(WARNING) << "Invalid number of elements in file: " << line; + } + continue; + } + + Cost value = .0; + if (!SimpleAtoValue(line, &value)) { + LOG(WARNING) << "Invalid value in file: " << line; + continue; + } + if (model.bin_capacity() <= .0) { + DCHECK_GT(value, .0); + model.set_bin_capacity(value); + } else { + model.AddItem(value); + } + } + DCHECK_GT(model.bin_capacity(), .0); + DCHECK_GT(model.weights().size(), .0); + DCHECK_EQ(num_items, model.weights().size()); + model.SortWeights(); + return model; +} + +BinPackingModel ReadCsp(absl::string_view filename) { + BinPackingModel model; + BaseInt num_item_types = 0; + for (const std::string& line : + FileLines(filename, FileLineIterator::REMOVE_INLINE_CR | + FileLineIterator::REMOVE_BLANK_LINES)) { + if (num_item_types == 0) { + if (!SimpleAtoValue(line, &num_item_types)) { + LOG(WARNING) << "Invalid number of elements in file: " << line; + } + continue; + } + if (model.bin_capacity() <= .0) { + Cost capacity = .0; + if (!SimpleAtoValue(line, &capacity)) { + LOG(WARNING) << "Invalid value in file: " << line; + } + model.set_bin_capacity(capacity); + continue; + } + + std::pair weight_and_demand = + absl::StrSplit(line, absl::ByAnyChar(" :\t"), absl::SkipEmpty()); + + Cost weight = .0; + if (!SimpleAtoValue(weight_and_demand.first, &weight)) { + LOG(WARNING) << "Invalid weight in file: " << line; + continue; + } + + BaseInt demand = 0; + if (!SimpleAtoValue(weight_and_demand.second, &demand)) { + LOG(WARNING) << "Invalid demand in file: " << line; + continue; + } + for (BaseInt i = 0; i < demand; ++i) { + model.AddItem(weight); + } + } + DCHECK_GT(model.bin_capacity(), .0); + DCHECK_GT(model.weights().size(), .0); + DCHECK_GE(num_item_types, model.num_items()); + model.SortWeights(); + return model; +} + +void BestFit(const ElementCostVector& weights, Cost bin_capacity, + const std::vector& items, PartialBins& bins_data) { + for (ElementIndex item : items) { + Cost item_weight = weights[item]; + BaseInt selected_bin = bins_data.bins.size(); + for (BaseInt bin = 0; bin < bins_data.bins.size(); ++bin) { + Cost max_load = bin_capacity - item_weight; + if (bins_data.loads[bin] <= max_load && + (selected_bin == bins_data.bins.size() || + bins_data.loads[bin] > bins_data.loads[selected_bin])) { + selected_bin = bin; + } + } + if (selected_bin == bins_data.bins.size()) { + bins_data.bins.emplace_back(); + bins_data.loads.emplace_back(); + } + bins_data.bins[selected_bin].push_back(item); + bins_data.loads[selected_bin] += item_weight; + } +} + +const SparseColumn& BinPackingSetCoverModel::BinPackingModelGlobals::GetBin( + SubsetIndex j) const { + if (j < SubsetIndex(full_model.num_subsets())) { + return full_model.columns()[j]; + } + DCHECK(candidate_bin != nullptr); + return *candidate_bin; +} + +uint64_t BinPackingSetCoverModel::BinHash::operator()(SubsetIndex j) const { + DCHECK(globals != nullptr); + return absl::HashOf(globals->GetBin(j)); +} + +uint64_t BinPackingSetCoverModel::BinEq::operator()(SubsetIndex j1, + SubsetIndex j2) const { + DCHECK(globals != nullptr); + return globals->GetBin(j1) == globals->GetBin(j2); +} + +bool BinPackingSetCoverModel::AddBin(const SparseColumn& bin) { + if (TryInsertBin(bin)) { + globals_.full_model.AddEmptySubset(1.0); + for (ElementIndex i : bin) { + globals_.full_model.AddElementToLastSubset(i); + } + return true; + } + return false; +} + +bool BinPackingSetCoverModel::TryInsertBin(const SparseColumn& bin) { + DCHECK(absl::c_is_sorted(bin)); + DCHECK(absl::c_adjacent_find(bin) == bin.end()); + DCHECK(globals_.candidate_bin == nullptr); + globals_.candidate_bin = &bin; + + SubsetIndex candidate_j(globals_.full_model.num_subsets()); + bool inserted = bin_set_.insert(candidate_j).second; + + globals_.candidate_bin = nullptr; + return inserted; +} + +void InsertBinsIntoModel(PartialBins& bins_data, + BinPackingSetCoverModel& model) { + for (BaseInt i = 0; i < bins_data.bins.size(); ++i) { + if (bins_data.bins[i].size() > 0) { + absl::c_sort(bins_data.bins[i]); + model.AddBin(bins_data.bins[i]); + } + } +} + +void AddRandomizedBins(const BinPackingModel& model, BaseInt num_bins, + BinPackingSetCoverModel& scp_model, std::mt19937& rnd) { + PartialBins bins_data; + std::vector items(model.num_items()); + absl::c_iota(items, ElementIndex(0)); + + while (scp_model.full_model().num_subsets() < num_bins) { + // Generate bins all containing a specific item + for (ElementIndex n : model.ItemRange()) { + BaseInt unique_bin_num = scp_model.full_model().num_subsets(); + VLOG_EVERY_N_SEC(1, 1) + << "[RGEN] Generating bins: " << unique_bin_num << " / " << num_bins + << " (" << 100.0 * unique_bin_num / num_bins << "%)"; + if (scp_model.full_model().num_subsets() >= num_bins) { + break; + } + + absl::c_shuffle(items, rnd); + + auto n_it = absl::c_find(items, n); + std::iter_swap(n_it, items.end() - 1); + items.pop_back(); + + bins_data.bins.clear(); + bins_data.loads.clear(); + for (BaseInt j = 0; j < 10; ++j) { + bins_data.bins.push_back({n}); + bins_data.loads.push_back(model.weights()[n]); + } + BestFit(model.weights(), model.bin_capacity(), items, bins_data); + InsertBinsIntoModel(bins_data, scp_model); + + items.push_back(n); + + if (unique_bin_num == scp_model.full_model().num_subsets()) { + VLOG(1) << "No new bins generated."; + break; + } + } + } + + scp_model.CompleteModel(); +} + +BinPackingSetCoverModel GenerateInitialBins(const BinPackingModel& model) { + BinPackingSetCoverModel scp_model(&model); + PartialBins bins_data; + std::vector items(model.num_items()); + + absl::c_iota(items, ElementIndex(0)); + absl::c_sort(items, [&](auto i1, auto i2) { + return model.weights()[i1] > model.weights()[i2]; + }); + BestFit(model.weights(), model.bin_capacity(), items, bins_data); + InsertBinsIntoModel(bins_data, scp_model); + VLOG(1) << "[BFIT] Largest first best-fit solution: " << bins_data.bins.size() + << " bins"; + + scp_model.CompleteModel(); + return scp_model; +} + +void ExpKnap::SaveBin() { + collected_bins_.back().clear(); + for (ElementIndex i : exceptions_) { + break_selection_[i] = !break_selection_[i]; + } + for (ElementIndex i : break_solution_) { + if (break_selection_[i]) { + collected_bins_.back().push_back(i); + } + } + for (ElementIndex i : exceptions_) { + if (break_selection_[i]) { + collected_bins_.back().push_back(i); + } + break_selection_[i] = !break_selection_[i]; + } + absl::c_sort(collected_bins_.back()); + DCHECK(absl::c_adjacent_find(collected_bins_.back()) == + collected_bins_.back().end()); +} +void ExpKnap::InitSolver(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit) { + capacity_ = capacity; + items_.resize(profits.size()); + collected_bins_.clear(); + for (ElementIndex i; i < ElementIndex(items_.size()); ++i) { + items_[i] = {profits[i], weights[i], i}; + } + absl::c_sort(items_, [](Item i1, Item i2) { + return i1.profit / i1.weight < i2.profit / i2.weight; + }); + + bnb_node_countdown_ = bnb_nodes_limit; + best_delta_ = .0; + exceptions_.clear(); + break_selection_.assign(profits.size(), false); + break_solution_.clear(); +} + +void ExpKnap::FindGoodColumns(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit) { + InitSolver(profits, weights, capacity, bnb_nodes_limit); + Cost curr_best_cost = .0; + PartialBins more_bins; + std::vector remaining_items; + + bnb_node_countdown_ = bnb_nodes_limit; + inserted_items_.assign(profits.size(), false); + do { + collected_bins_.emplace_back(); + Heuristic(); + VLOG(5) << "[KPCG] Heuristic solution: cost " + << break_profit_sum_ + best_delta_; + EleBranch(); + + for (ElementIndex i : collected_bins_.back()) { + inserted_items_[i] = true; + } + gtl::STLEraseAllFromSequenceIf( + &items_, [&](Item item) { return inserted_items_[item.index]; }); + } while (!items_.empty() && break_profit_sum_ + best_delta_ > 1.0); +} + +bool ExpKnap::EleBranch() { + exceptions_.clear(); + return EleBranch(.0, break_weight_sum_ - capacity_, break_it_ - 1, break_it_); +} + +namespace { +static Cost BoundCheck(Cost best_delta, Cost profit_delta, Cost overweight, + ExpKnap::Item item) { + Cost bound = best_delta + 0.0; // + 1.0 for integral profits + return (profit_delta - bound) * item.weight - overweight * item.profit; +} +} // namespace + +// Adapted version from David Pisinger elebranch function: +// https://hjemmesider.diku.dk/~pisinger/expknap.c +bool ExpKnap::EleBranch(Cost profit_delta, Cost overweigth, ItemIt out_item, + ItemIt in_item) { + VLOG(6) << "[KPCG] EleBranch: profit_delta " << profit_delta << " overweigth " + << overweigth; + if (bnb_node_countdown_-- <= 0) { + return false; + } + bool improved = false; + + if (overweigth <= .0) { + if (profit_delta > best_delta_) { + best_delta_ = profit_delta; + improved = true; + VLOG(5) << "[KPCG] Improved best cost " + << break_profit_sum_ + best_delta_; + } + + bool maximal = true; + while (bnb_node_countdown_ > 0 && in_item < items_.rend() && + BoundCheck(best_delta_, profit_delta, overweigth, *in_item) >= 0) { + exceptions_.push_back(in_item->index); + Cost next_delta = profit_delta + in_item->profit; + Cost next_oweight = overweigth + in_item->weight; + maximal &= !EleBranch(next_delta, next_oweight, out_item, ++in_item); + exceptions_.pop_back(); + } + + if (improved && maximal) { + SaveBin(); + } + improved |= !maximal; + } else { + while (bnb_node_countdown_ > 0 && out_item >= items_.rbegin() && + BoundCheck(best_delta_, profit_delta, overweigth, *out_item) >= 0) { + exceptions_.push_back(out_item->index); + Cost next_delta = profit_delta - out_item->profit; + Cost next_oweight = overweigth - out_item->weight; + improved |= EleBranch(next_delta, next_oweight, --out_item, in_item); + exceptions_.pop_back(); + } + } + return improved; +} + +void ExpKnap::Heuristic() { + best_delta_ = break_profit_sum_ = break_weight_sum_ = .0; + break_it_ = items_.rbegin(); + break_selection_.assign(items_.size(), false); + break_solution_.clear(); + exceptions_.clear(); + while (break_it_ < items_.rend() && + break_it_->weight <= capacity_ - break_weight_sum_) { + break_profit_sum_ += break_it_->profit; + break_weight_sum_ += break_it_->weight; + break_selection_[break_it_->index] = true; + break_solution_.push_back(break_it_->index); + ++break_it_; + } + Cost residual = capacity_ - break_weight_sum_; + + VLOG(5) << "[KPCG] Break solution: cost " << break_profit_sum_ + << ", residual " << residual; + + Cost profit_delta_ub = residual * break_it_->profit / break_it_->weight; + if (profit_delta_ub == .0) { + SaveBin(); + return; + } + + // Try filling the residual space with less efficient (maybe smaller) items + best_delta_ = .0; + for (auto it = break_it_; it < items_.rend(); it++) { + if (it->weight <= residual && it->profit > best_delta_) { + exceptions_ = {it->index}; + best_delta_ = it->profit; + if (best_delta_ >= profit_delta_ub) { + SaveBin(); + return; + } + } + } + + // Try removing an item and adding the break item + Cost min_weight = break_it_->weight - residual; + for (auto it = break_it_ - 1; it >= items_.rbegin(); it--) { + Cost profit_delta = break_it_->profit - it->profit; + if (it->weight >= min_weight && profit_delta > best_delta_) { + exceptions_ = {break_it_->index, it->index}; + best_delta_ = profit_delta; + if (best_delta_ >= profit_delta_ub) { + SaveBin(); + return; + } + } + } + SaveBin(); +} + +bool BinPackingSetCoverModel::UpdateCore( + Cost best_lower_bound, const ElementCostVector& best_multipliers, + const scp::Solution& best_solution, bool force) { + if (!base::IsTimeToUpdate(best_lower_bound, force)) { + return false; + } + + Cost full_lower_bound = base::UpdateMultipliers(best_multipliers); + if (scp::DivideIfGE0(std::abs(full_lower_bound - best_lower_bound), + best_lower_bound) < 0.01 && + best_lower_bound != prev_lower_bound_) { + prev_lower_bound_ = best_lower_bound; + knapsack_solver_.FindGoodColumns(best_multipliers, bpp_model_->weights(), + bpp_model_->bin_capacity(), + /*bnb_nodes_limit=*/1000); + BaseInt num_added_bins = 0; + for (const SparseColumn& bin : knapsack_solver_.collected_bins()) { + num_added_bins += AddBin(bin) ? 1 : 0; + } + if (num_added_bins > 0) { + VLOG(4) << "[KPCG] Added " << num_added_bins << " / " + << globals_.full_model.num_subsets() << " bins"; + } + base::SizeUpdate(); + // TODO(user): add incremental update only for the new columns just added + base::UpdateMultipliers(best_multipliers); + } + + if (base::num_focus_subsets() < FixingFullModelView().num_focus_subsets()) { + ComputeAndSetFocus(best_lower_bound, best_solution); + return true; + } + return false; +} +} // namespace operations_research diff --git a/ortools/algorithms/bin_packing.h b/ortools/algorithms/bin_packing.h new file mode 100644 index 00000000000..17ce1417841 --- /dev/null +++ b/ortools/algorithms/bin_packing.h @@ -0,0 +1,180 @@ + +// Copyright 2025 Francesco Cavaliere +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H +#define OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H + +#include +#include +#include +#include +#include + +#include +#include + +#include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_cft.h" + +namespace operations_research { + +class BinPackingModel { + public: + BinPackingModel() = default; + BaseInt num_items() const { return weigths_.size(); } + Cost bin_capacity() const { return bin_capcaity_; } + void set_bin_capacity(Cost capacity); + const ElementCostVector& weights() const { return weigths_; } + void AddItem(Cost weight); + void SortWeights(); + ElementRange ItemRange() const { + return {ElementIndex(), ElementIndex(weigths_.size())}; + } + + private: + bool is_sorted_ = false; + Cost bin_capcaity_ = .0; + ElementCostVector weigths_ = {}; +}; + +struct PartialBins { + std::vector bins; + std::vector loads; +}; + +using SubsetHashVector = util_intops::StrongVector; + +class ExpKnap { + public: + struct Item { + Cost profit; // profit + Cost weight; // weight + ElementIndex index; + }; + using ItemIt = std::vector::const_reverse_iterator; + + void SaveBin(); + void FindGoodColumns(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit); + void InitSolver(const ElementCostVector& profits, + const ElementCostVector& weights, Cost capacity, + BaseInt bnb_nodes_limit); + + bool EleBranch(); + bool EleBranch(Cost profit_sum, Cost overweight, ItemIt out_item, + ItemIt in_item); + + void Heuristic(); + + const std::vector& collected_bins() const { + return collected_bins_; + } + Cost best_cost() const { return break_profit_sum_ + best_delta_; } + + private: + Cost capacity_; // capacity + util_intops::StrongVector items_; // items + ItemIt break_it_; + Cost break_profit_sum_; + Cost break_weight_sum_; + Cost best_delta_; + std::vector exceptions_; + std::vector collected_bins_; + ElementBoolVector break_selection_; + ElementBoolVector inserted_items_; + SparseColumn break_solution_; + BaseInt bnb_node_countdown_; + std::mt19937 rnd_; +}; + +class BinPackingSetCoverModel : public scp::FullToCoreModel { + using base = scp::FullToCoreModel; + + struct BinPackingModelGlobals { + // Dirty hack to avoid invalidation of pointers/references + // A pointer to this data structure is used to compute the hash of bins + // starting from their indices. + scp::Model full_model; + + // External bins do not have a valid index in the model, a temporary pointer + // is used insteda (even dirtier hack). + const SparseColumn* candidate_bin; + + const SparseColumn& GetBin(SubsetIndex j) const; + }; + + struct BinHash { + const BinPackingModelGlobals* globals; + uint64_t operator()(SubsetIndex j) const; + }; + + struct BinEq { + const BinPackingModelGlobals* globals; + uint64_t operator()(SubsetIndex j1, SubsetIndex j2) const; + }; + + public: + BinPackingSetCoverModel(const BinPackingModel* bpp_model) + : globals_{scp::Model(), nullptr}, + bpp_model_(bpp_model), + knapsack_solver_(), + bin_set_({}, 0, BinHash{&globals_}, BinEq{&globals_}), + column_gen_countdown_(10), + column_gen_period_(10) {} + const scp::Model& full_model() const { return globals_.full_model; } + + bool AddBin(const SparseColumn& bin); + void CompleteModel() { + globals_.full_model.CreateSparseRowView(); + static_cast(*this) = + scp::FullToCoreModel(&globals_.full_model); + } + + bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const scp::Solution& best_solution, bool force) override; + + private: + bool TryInsertBin(const SparseColumn& bin); + + BinPackingModelGlobals globals_; + const BinPackingModel* bpp_model_; + ExpKnap knapsack_solver_; + + // Contains bin indices, but it really should contains "bins" (aka, + // SparseColumn). However to avoid redundant allocations (in scp::Model and + // in the set) we cannot store them also here. We cannot also use + // iterators/pointers/references because they can be invalidated. So we store + // bin indices and do ungodly hacky shenanigans to get the bins from them. + absl::flat_hash_set bin_set_; + + Cost prev_lower_bound_; + BaseInt column_gen_countdown_; + BaseInt column_gen_period_; +}; + +BinPackingModel ReadBpp(absl::string_view filename); +BinPackingModel ReadCsp(absl::string_view filename); + +void BestFit(const BinPackingModel& model, + const std::vector& items, PartialBins& bins_data); + +BinPackingSetCoverModel GenerateInitialBins(const BinPackingModel& model); + +void AddRandomizedBins(const BinPackingModel& model, BaseInt num_bins, + BinPackingSetCoverModel& scp_model, std::mt19937& rnd); + +} // namespace operations_research +#endif /* OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H */ diff --git a/ortools/algorithms/samples/bin_packing_cft.cc b/ortools/algorithms/samples/bin_packing_cft.cc new file mode 100644 index 00000000000..0eafd474160 --- /dev/null +++ b/ortools/algorithms/samples/bin_packing_cft.cc @@ -0,0 +1,144 @@ + +// Copyright 2025 Francesco Cavaliere +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include +#include + +#include +#include + +#include "ortools/algorithms/bin_packing.h" +#include "ortools/base/init_google.h" +#include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_cft.h" + +using namespace operations_research; +ABSL_FLAG(std::string, instance, "", "BPP instance int RAIL format."); +ABSL_FLAG(int, bins, 1000, "Number of bins to generate."); + +template +std::string Stringify(const Iterable& col) { + std::string result; + for (auto i : col) { + absl::StrAppend(&result, " ", i); + } + return result; +} + +bool operator==(const SparseColumn& lhs, const std::vector& rhs) { + if (lhs.size() != rhs.size()) return false; + auto lit = lhs.begin(); + auto rit = rhs.begin(); + while (lit != lhs.end() && rit != rhs.end()) { + if (static_cast(*lit) != *rit) { + return false; + } + ++lit; + ++rit; + } + return true; +} + +void RunTest(const ElementCostVector& weights, const ElementCostVector& profits, + const std::vector& expected) { + ExpKnap knap_solver; + + for (ElementIndex i; i < ElementIndex(weights.size()); ++i) { + std::cout << "Item " << i << " -- profit: " << profits[i] + << " weight: " << weights[i] + << " efficiency: " << profits[i] / weights[i] << "\n"; + } + + knap_solver.InitSolver(profits, weights, 6, 100000000); + knap_solver.Heuristic(); + std::cout << "Heur solution cost " << knap_solver.best_cost() << " -- " + << Stringify(knap_solver.collected_bins()[0]) << "\n"; + + knap_solver.EleBranch(); + std::cout << "B&b solution cost " << knap_solver.best_cost() << " -- " + << Stringify(knap_solver.collected_bins()[0]) << "\n"; + + const auto& result = knap_solver.collected_bins()[0]; + if (!(result == expected)) { + std::cout << "Error: expected " << Stringify(expected) << " but got " + << Stringify(result) << "\n"; + } + std::cout << std::endl; +} + +void KnapsackTest() { + std::cout << "Testing knapsack\n"; + ExpKnap knap_solver; + ElementCostVector ws = {1, 2, 3, 4, 5}; + RunTest(ws, {10, 20, 30, 40, 51}, {0, 4}); + RunTest(ws, {10, 20, 30, 41, 50}, {1, 3}); + RunTest(ws, {10, 20, 31, 40, 50}, {0, 1, 2}); + RunTest(ws, {10, 21, 30, 41, 50}, {1, 3}); + RunTest(ws, {11, 21, 30, 40, 50}, {0, 1, 2}); + RunTest(ws, {11, 20, 31, 40, 50}, {0, 1, 2}); + RunTest(ws, {11, 20, 30, 41, 50}, {0, 4}); + RunTest(ws, {11, 20, 30, 40, 51}, {0, 4}); + RunTest(ws, {11, 21, 31, 40, 50}, {0, 1, 2}); + RunTest({4.1, 2, 2, 2}, {8.5, 3, 3, 3}, {1, 2, 3}); +} + +int main(int argc, char** argv) { + InitGoogle(argv[0], &argc, &argv, true); + + // KnapsackTest(); + // return 0; + + BinPackingModel model = ReadBpp(absl::GetFlag(FLAGS_instance)); + + // Quick run with a minimal set of bins + BinPackingSetCoverModel scp_model = GenerateInitialBins(model); + scp::PrimalDualState best_result = scp::RunCftHeuristic(scp_model); + + if (absl::GetFlag(FLAGS_bins) > 0) { + // Run the CFT again with more bins to get a better solution + std::mt19937 rnd(0); + AddRandomizedBins(model, absl::GetFlag(FLAGS_bins), scp_model, rnd); + scp::PrimalDualState result = + scp::RunCftHeuristic(scp_model, best_result.solution); + if (result.solution.cost() < best_result.solution.cost()) { + best_result = result; + } + } + + auto [solution, dual] = best_result; + if (solution.subsets().empty()) { + std::cerr << "Error: failed to find any solution\n"; + } else { + std::cout << "Solution: " << solution.cost() << '\n'; + } + + if (dual.multipliers().empty()) { + std::cerr << "Error: failed to find any dual\n"; + } else { + std::cout << "Core Lower bound: " << dual.lower_bound() << '\n'; + } + + // The lower bound computed on the full model is not a real lower bound unless + // the knapsack subproblem failed to fined any negative reduced cost bin to + // add to the set cover model. + // TODO(anyone): add a flag to indicate if a valid LB has been found or not. + if (scp_model.best_dual_state().multipliers().empty()) { + std::cerr << "Error: no real dual state has been computed\n"; + } else { + std::cout << "Restricted Lower bound: " + << scp_model.best_dual_state().lower_bound() << '\n'; + } + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/ortools/set_cover/samples/cft.cc b/ortools/set_cover/samples/cft.cc new file mode 100644 index 00000000000..90a72ea6b0b --- /dev/null +++ b/ortools/set_cover/samples/cft.cc @@ -0,0 +1,53 @@ +#include +#include +#include + +#include "ortools/base/init_google.h" +#include "ortools/set_cover/set_cover_cft.h" +#include "ortools/set_cover/set_cover_reader.h" + +using namespace operations_research; +ABSL_FLAG(std::string, instance, "", "SCP instance int RAIL format."); + +#define DO_PRICING +int main(int argc, char **argv) { + InitGoogle(argv[0], &argc, &argv, true); + + scp::Model original_model = ReadOrlibRail(absl::GetFlag(FLAGS_instance)); + +#ifdef DO_PRICING + scp::FullToCoreModel model(&original_model); +#else + scp::SubModel model(&original_model); +#endif + + scp::PrimalDualState result = scp::RunCftHeuristic(model); + auto [solution, dual] = result; + if (solution.subsets().empty()) { + std::cerr << "Error: failed to find any solution\n"; + } else { + std::cout << "Solution: " << solution.cost() << '\n'; + } + +#ifdef DO_PRICING + if (dual.multipliers().empty()) { + std::cerr << "Error: failed to find any dual\n"; + } else { + std::cout << "Core Lower bound: " << dual.lower_bound() << '\n'; + } + if (model.best_dual_state().multipliers().empty()) { + std::cerr << "Error: no real dual state has been computed\n"; + } else { + std::cout << "Full Lower bound: " << model.best_dual_state().lower_bound() + << '\n'; + } +#else + if (dual.multipliers().empty()) { + std::cerr << "Error: failed to find any dual\n"; + } else { + std::cout << "Lower bound: " << dual.lower_bound() << '\n'; + } +#endif + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/ortools/set_cover/set_cover_cft.cc b/ortools/set_cover/set_cover_cft.cc index 71628ed1b30..88af8d42d57 100644 --- a/ortools/set_cover/set_cover_cft.cc +++ b/ortools/set_cover/set_cover_cft.cc @@ -14,39 +14,73 @@ #include "ortools/set_cover/set_cover_cft.h" #include +#include +#include #include #include #include +#include #include +#include #include #include -#include "ortools/base/status_macros.h" #include "ortools/base/stl_util.h" #include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_submodel.h" #include "ortools/set_cover/set_cover_views.h" +#include "ortools/set_cover/views.h" namespace operations_research::scp { +// Minimum distance between lower and upper bounds to consider them different. +// If cost are all integral, can be set neear to 1.0 +#define CFT_BOUND_EPSILON .999 +#define CFT_MAX_MULTIPLIER 1e9 +#define CFT_MEASURE_TIME + //////////////////////////////////////////////////////////////////////// ////////////////////////// COMMON DEFINITIONS ////////////////////////// //////////////////////////////////////////////////////////////////////// namespace { -class CoverCounters { + +#ifdef CFT_MEASURE_TIME +#define CFT_MEASURE_SCOPE_DURATION(Timer) \ + Timer.Start(); \ + Defer pause_timer = [&] { Timer.Stop(); }; + +thread_local WallTimer subgradient_time; +thread_local WallTimer greedy_time; +thread_local WallTimer three_phase_time; +thread_local WallTimer refinement_time; +#else +#define CFT_MEASURE_SCOPE_DURATION(Timer) +#endif + +// StopWatch does not add up duration of multiple invocations, Defer is a lower +// level construct useful in this case. +template +struct Defer : T { + Defer(T&& t) : T(std::forward(t)) {} + ~Defer() { T::operator()(); } +}; + +template +class CoverCountersImpl { public: - CoverCounters(BaseInt nelems = 0) : cov_counters(nelems, 0) {} + CoverCountersImpl(BaseInt nelems = 0) : cov_counters(nelems, 0) {} void Reset(BaseInt nelems) { cov_counters.assign(nelems, 0); } BaseInt Size() const { return cov_counters.size(); } - BaseInt operator[](ElementIndex i) const { - assert(i < ElementIndex(cov_counters.size())); + BaseInt operator[](IndexT i) const { + assert(i < IndexT(cov_counters.size())); return cov_counters[i]; } template BaseInt Cover(const IterableT& subset) { BaseInt covered = 0; - for (ElementIndex i : subset) { + for (IndexT i : subset) { covered += cov_counters[i] == 0 ? 1ULL : 0ULL; cov_counters[i]++; } @@ -56,7 +90,7 @@ class CoverCounters { template BaseInt Uncover(const IterableT& subset) { BaseInt uncovered = 0; - for (ElementIndex i : subset) { + for (IndexT i : subset) { --cov_counters[i]; uncovered += cov_counters[i] == 0 ? 1ULL : 0ULL; } @@ -67,19 +101,22 @@ class CoverCounters { template bool IsRedundantCover(IterableT const& subset) const { return absl::c_all_of(subset, - [&](ElementIndex i) { return cov_counters[i] > 0; }); + [&](IndexT i) { return cov_counters[i] > 0; }); } // Check if all the elements would still be covered if the subset was removed. template bool IsRedundantUncover(IterableT const& subset) const { return absl::c_all_of(subset, - [&](ElementIndex i) { return cov_counters[i] > 1; }); + [&](IndexT i) { return cov_counters[i] > 1; }); } private: - ElementToIntVector cov_counters; + util_intops::StrongVector cov_counters; }; + +using CoverCounters = CoverCountersImpl; +using FullCoverCounters = CoverCountersImpl; } // namespace Solution::Solution(const SubModel& model, @@ -106,9 +143,10 @@ BoundCBs::BoundCBs(const SubModel& model) prev_best_lb_(std::numeric_limits::lowest()), max_iter_countdown_(10 * model.num_focus_elements()), // Arbitrary from [1] - exit_test_countdown_(300), // Arbitrrary from [1] + exit_test_countdown_(300), // Arbitrary from [1] exit_test_period_(300), // Arbitrary from [1] - step_size_(0.1), // Arbitrary from [1] + unfixed_run_extension_(0), + step_size_(0.1), // Arbitrary from [1] last_min_lb_seen_(std::numeric_limits::max()), last_max_lb_seen_(.0), step_size_update_countdown_(20), // Arbitrary from [1] @@ -116,61 +154,71 @@ BoundCBs::BoundCBs(const SubModel& model) {} bool BoundCBs::ExitCondition(const SubgradientContext& context) { - if (last_core_update_countdown_ > 0) { - --last_core_update_countdown_; - return false; - } - if (--max_iter_countdown_ <= 0 || squared_norm_ <= kTol || - // Note: assumes integral costs - context.best_dual_state.lower_bound() >= - context.best_solution.cost() - context.model.fixed_cost() - .999) { + Cost best_lb = context.best_lower_bound; + Cost best_ub = context.best_solution.cost() - context.model.fixed_cost(); + if (--max_iter_countdown_ <= 0 || squared_norm_ <= kTol) { return true; } if (--exit_test_countdown_ > 0) { return false; } + if (prev_best_lb_ >= best_ub - CFT_BOUND_EPSILON) { + return true; + } exit_test_countdown_ = exit_test_period_; - Cost best_lower_bound = context.best_dual_state.lower_bound(); - Cost abs_improvement = best_lower_bound - prev_best_lb_; - Cost rel_improvement = DivideIfGE0(abs_improvement, best_lower_bound); - prev_best_lb_ = best_lower_bound; - return abs_improvement < 1.0 && rel_improvement < .001; + Cost abs_improvement = best_lb - prev_best_lb_; + Cost rel_improvement = DivideIfGE0(abs_improvement, best_lb); + prev_best_lb_ = best_lb; + + if (abs_improvement >= 1.0 || rel_improvement >= .001) { + return false; + } + + // (Not in [1]): During the first unfixed iteration we want to converge closer + // to the optimum + VLOG(3) << "[SUBG] First iteration extension " << unfixed_run_extension_; + BaseInt extension = (context.model.fixed_columns().empty() ? 2 : 1); + return unfixed_run_extension_++ >= extension; } void BoundCBs::ComputeMultipliersDelta(const SubgradientContext& context, ElementCostVector& delta_mults) { - Cost lower_bound = context.current_dual_state.lower_bound(); - UpdateStepSize(lower_bound); - direction_ = context.subgradient; - MakeMinimalCoverageSubgradient(context, direction_); - - if (prev_direction_.empty()) { - prev_direction_ = direction_; + BaseInt extension = (context.model.fixed_columns().empty() ? 2 : 1); + if (unfixed_run_extension_ < extension) { + MakeMinimalCoverageSubgradient(context, direction_); } squared_norm_ = .0; - // Stabilize current direction and compute squared norm for (ElementIndex i : context.model.ElementRange()) { - direction_[i] = stabilization_coeff * direction_[i] + - (1.0 - stabilization_coeff) * prev_direction_[i]; - if ((context.current_dual_state.multipliers()[i] <= .0) && - (direction_[i] < .0)) { + Cost curr_i_mult = context.current_dual_state.multipliers()[i]; + if ((curr_i_mult <= .0 && direction_[i] < .0) || + (curr_i_mult > CFT_MAX_MULTIPLIER && direction_[i] > .0)) { direction_[i] = 0; } + squared_norm_ += direction_[i] * direction_[i]; - prev_direction_[i] = direction_[i]; } + if (squared_norm_ <= kTol) { + delta_mults.assign(context.model.num_elements(), .0); + return; + } + + UpdateStepSize(context); Cost upper_bound = context.best_solution.cost() - context.model.fixed_cost(); + Cost lower_bound = context.current_dual_state.lower_bound(); Cost delta = upper_bound - lower_bound; - Cost step_constant = step_size_ * delta / squared_norm_; + Cost step_constant = (step_size_ * delta) / squared_norm_; + for (ElementIndex i : context.model.ElementRange()) { delta_mults[i] = step_constant * direction_[i]; + DCHECK(std::isfinite(delta_mults[i])); } } -void BoundCBs::UpdateStepSize(Cost lower_bound) { +void BoundCBs::UpdateStepSize(SubgradientContext context) { + Cost lower_bound = context.current_dual_state.lower_bound(); last_min_lb_seen_ = std::min(last_min_lb_seen_, lower_bound); last_max_lb_seen_ = std::max(last_max_lb_seen_, lower_bound); @@ -181,12 +229,10 @@ void BoundCBs::UpdateStepSize(Cost lower_bound) { Cost gap = DivideIfGE0(delta, last_max_lb_seen_); if (gap <= .001) { // Arbitray from [1] step_size_ *= 1.5; // Arbitray from [1] - - // Arbitrary from c4v4 - stabilization_coeff = (1.0 + stabilization_coeff) / 2.0; - + VLOG(4) << "[SUBG] Sep size set at " << step_size_; } else if (gap > .01) { // Arbitray from [1] step_size_ /= 2.0; // Arbitray from [1] + VLOG(4) << "[SUBG] Sep size set at " << step_size_; } last_min_lb_seen_ = std::numeric_limits::max(); last_max_lb_seen_ = .0; @@ -217,15 +263,15 @@ void BoundCBs::MakeMinimalCoverageSubgradient(const SubgradientContext& context, } } -bool BoundCBs::UpdateCoreModel(SubModel& core_model, - PrimalDualState& best_state) { - if (core_model.UpdateCore(best_state)) { - // grant at least 10 iterations before the next exit test - prev_best_lb_ = - std::min(prev_best_lb_, best_state.dual_state.lower_bound()); - exit_test_countdown_ = std::max(exit_test_countdown_, 20); - max_iter_countdown_ = std::max(max_iter_countdown_, 20); - last_core_update_countdown_ = 20; +bool BoundCBs::UpdateCoreModel(SubgradientContext context, + CoreModel& core_model, bool force) { + if (core_model.UpdateCore(context.best_lower_bound, context.best_multipliers, + context.best_solution, force)) { + prev_best_lb_ = std::numeric_limits::lowest(); + // Grant at least `min_iters` iterations before the next exit test + constexpr BaseInt min_iters = 10; + exit_test_countdown_ = std::max(exit_test_countdown_, min_iters); + max_iter_countdown_ = std::max(max_iter_countdown_, min_iters); return true; } return false; @@ -233,20 +279,33 @@ bool BoundCBs::UpdateCoreModel(SubModel& core_model, void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, PrimalDualState& best_state) { + CFT_MEASURE_SCOPE_DURATION(subgradient_time); DCHECK(ValidateSubModel(model)); ElementCostVector subgradient = ElementCostVector(model.num_elements(), .0); DualState dual_state = best_state.dual_state; + Cost best_lower_bound = dual_state.lower_bound(); + ElementCostVector best_multipliers = dual_state.multipliers(); Solution solution; ElementCostVector multipliers_delta(model.num_elements()); // to avoid allocs SubgradientContext context = {.model = model, .current_dual_state = dual_state, - .best_dual_state = best_state.dual_state, + .best_lower_bound = best_lower_bound, + .best_multipliers = best_multipliers, .best_solution = best_state.solution, .subgradient = subgradient}; - size_t iter = 0; - while (!cbs.ExitCondition(context)) { + + for (size_t iter = 1; !cbs.ExitCondition(context); ++iter) { + // Poor multipliers can lead to wasted iterations or stagnation in the + // subgradient method. To address this, we adjust the multipliers to + // get closer the trivial lower bound (= 0). + if (dual_state.lower_bound() < .0) { + VLOG(4) << "[SUBG] Dividing multipliers by 10"; + dual_state.DualUpdate( + model, [&](ElementIndex i, Cost& i_mult) { i_mult /= 10.0; }); + } + // Compute subgradient (O(nnz)) subgradient.assign(model.num_elements(), 1.0); for (SubsetIndex j : model.SubsetRange()) { @@ -259,10 +318,12 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, cbs.ComputeMultipliersDelta(context, multipliers_delta); dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { - i_mult = std::max(.0, i_mult + multipliers_delta[i]); + i_mult = std::clamp(i_mult + multipliers_delta[i], .0, + CFT_MAX_MULTIPLIER); }); - if (dual_state.lower_bound() > best_state.dual_state.lower_bound()) { - best_state.dual_state = dual_state; + if (dual_state.lower_bound() > best_lower_bound) { + best_lower_bound = dual_state.lower_bound(); + best_multipliers = dual_state.multipliers(); } cbs.RunHeuristic(context, solution); @@ -271,23 +332,35 @@ void SubgradientOptimization(SubModel& model, SubgradientCBs& cbs, best_state.solution = solution; } - if (++iter % 50 == 0) - std::cout << "Subgradient " << iter << " -- Bounds: Lower " - << dual_state.lower_bound() << ", best " - << best_state.dual_state.lower_bound() << " - Upper " - << best_state.solution.cost() - model.fixed_cost() - << ", global " << best_state.solution.cost() << "\n"; + VLOG_EVERY_N(4, 100) << "[SUBG] " << iter << ": Bounds: Lower " + << dual_state.lower_bound() << ", best " + << best_lower_bound << " - Upper " + << best_state.solution.cost() - model.fixed_cost() + << ", global " << best_state.solution.cost(); - if (cbs.UpdateCoreModel(model, best_state)) { - std::cout << "Core model has been updated.\n"; - dual_state = best_state.dual_state; + if (cbs.UpdateCoreModel(context, model)) { + dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { + i_mult = best_multipliers[i]; + }); + best_lower_bound = dual_state.lower_bound(); } } - std::cout << "Subgradient End" << " -- Bounds: Lower " - << dual_state.lower_bound() << ", best " - << best_state.dual_state.lower_bound() << " - Upper " - << best_state.solution.cost() - model.fixed_cost() << ", global " - << best_state.solution.cost() << "\n"; + + if (cbs.UpdateCoreModel(context, model, /*force=*/true)) { + dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { + i_mult = best_multipliers[i]; + }); + best_lower_bound = dual_state.lower_bound(); + } + best_state.dual_state.DualUpdate(model, [&](ElementIndex i, Cost& i_mult) { + i_mult = best_multipliers[i]; + }); + DCHECK_EQ(best_state.dual_state.lower_bound(), best_lower_bound); + + VLOG(3) << "[SUBG] End - Bounds: Lower " << dual_state.lower_bound() + << ", best " << best_lower_bound << " - Upper " + << best_state.solution.cost() - model.fixed_cost() << ", global " + << best_state.solution.cost(); } //////////////////////////////////////////////////////////////////////// @@ -544,6 +617,8 @@ Solution RunMultiplierBasedGreedy(const SubModel& model, Cost CoverGreedly(const SubModel& model, const DualState& dual_state, Cost cost_cutoff, BaseInt stop_size, std::vector& sol_subsets) { + CFT_MEASURE_SCOPE_DURATION(greedy_time); + Cost sol_cost = .0; for (SubsetIndex j : sol_subsets) { sol_cost += model.subset_costs()[j]; @@ -588,6 +663,7 @@ Cost CoverGreedly(const SubModel& model, const DualState& dual_state, // Either remove redundant columns or discard solution RedundancyRemover remover(model, covered_rows); // TODO(?): cache it! + return remover.TryRemoveRedundantCols(model, cost_cutoff, sol_subsets); } @@ -648,12 +724,12 @@ void FixBestColumns(SubModel& model, PrimalDualState& state) { fix_at_least, cols_to_fix); // Fix columns and update the model - Cost fixed_cost_delta = model.FixColumns(cols_to_fix); + Cost fixed_cost_delta = model.FixMoreColumns(cols_to_fix); - std::cout << "Fixed " << cols_to_fix.size() - << " new columns with cost: " << fixed_cost_delta << '\n'; - std::cout << "Globally fixed " << model.fixed_columns().size() - << " columns, with cost " << model.fixed_cost() << '\n'; + VLOG(3) << "[3FIX] Fixed " << cols_to_fix.size() + << " new columns with cost: " << fixed_cost_delta; + VLOG(3) << "[3FIX] Globally fixed " << model.fixed_columns().size() + << " columns, with cost " << model.fixed_cost(); // Update multipliers for the reduced model dual_state.DualUpdate(model, [&](ElementIndex core_i, Cost& multiplier) { @@ -665,6 +741,9 @@ void FixBestColumns(SubModel& model, PrimalDualState& state) { void RandomizeDualState(const SubModel& model, DualState& dual_state, std::mt19937& rnd) { + // In [1] this step is described, not completely sure if it actually helps + // or not. Seems to me one of those "throw in some randomness, it never hurts" + // thing. dual_state.DualUpdate(model, [&](ElementIndex, Cost& i_multiplier) { i_multiplier *= absl::Uniform(rnd, 0.9, 1.1); }); @@ -696,6 +775,7 @@ void HeuristicCBs::ComputeMultipliersDelta(const SubgradientContext& context, } PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { + CFT_MEASURE_SCOPE_DURATION(three_phase_time); DCHECK(ValidateSubModel(model)); PrimalDualState best_state = {.solution = init_solution, @@ -704,52 +784,188 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { best_state.solution = RunMultiplierBasedGreedy(model, best_state.dual_state); } - std::cout << "Initial lower bound: " << best_state.dual_state.lower_bound() - << "\nInitial solution cost: " << best_state.solution.cost() - << "\nStarting 3-phase algorithm\n"; + VLOG(2) << "[3PHS] Initial lower bound: " + << best_state.dual_state.lower_bound() + << ", Initial solution cost: " << best_state.solution.cost() + << ", Starting 3-phase algorithm"; PrimalDualState curr_state = best_state; BaseInt iter_count = 0; std::mt19937 rnd(0xcf7); - while (model.num_focus_elements() > 0 && - // note: assumes integral costs - curr_state.dual_state.lower_bound() < - best_state.solution.cost() - model.fixed_cost() - .999) { + while (model.num_focus_elements() > 0) { ++iter_count; - std::cout << "\n\n3Phase iteration: " << iter_count << '\n'; - std::cout << " Active size: rows " << model.num_focus_elements() << "/" - << model.num_elements() << ", columns " - << model.num_focus_subsets() << "/" << model.num_subsets() - << '\n'; + VLOG(2) << "[3PHS] 3Phase iteration: " << iter_count; + VLOG(2) << "[3PHS] Active size: rows " << model.num_focus_elements() << "/" + << model.num_elements() << ", columns " << model.num_focus_subsets() + << "/" << model.num_subsets(); // Phase 1: refine the current dual_state and model BoundCBs dual_bound_cbs(model); - std::cout << "\nSubgradient Phase:\n"; + VLOG(2) << "[3PHS] Subgradient Phase:"; SubgradientOptimization(model, dual_bound_cbs, curr_state); if (iter_count == 1) { best_state.dual_state = curr_state.dual_state; } + if (curr_state.dual_state.lower_bound() >= + best_state.solution.cost() - model.fixed_cost() - CFT_BOUND_EPSILON) { + break; + } + // Phase 2: search for good solutions HeuristicCBs heuristic_cbs; heuristic_cbs.set_step_size(dual_bound_cbs.step_size()); - std::cout << "\nHeuristic Phase:\n"; + VLOG(2) << "[3PHS] Heuristic Phase:"; SubgradientOptimization(model, heuristic_cbs, curr_state); + if (iter_count == 1 && best_state.dual_state.lower_bound() < + curr_state.dual_state.lower_bound()) { + best_state.dual_state = curr_state.dual_state; + } if (curr_state.solution.cost() < best_state.solution.cost()) { best_state.solution = curr_state.solution; } + if (curr_state.dual_state.lower_bound() >= + best_state.solution.cost() - model.fixed_cost() - CFT_BOUND_EPSILON) { + break; + } - std::cout << "\n3Phase Bounds: Lower " - << best_state.dual_state.lower_bound() << ", Upper " - << best_state.solution.cost() << '\n'; + VLOG(2) << "[3PHS] 3Phase Bounds: Residual Lower " + << curr_state.dual_state.lower_bound() + model.fixed_cost() + << ", Upper " << best_state.solution.cost(); // Phase 3: Fix the best columns (diving) FixBestColumns(model, curr_state); RandomizeDualState(model, curr_state.dual_state, rnd); } - std::cout << "\n\n3Phase End\n"; - std::cout << "Final Bounds: Lower " << best_state.dual_state.lower_bound() - << ", Upper " << best_state.solution.cost() << '\n'; + VLOG(2) << "[3PHS] 3Phase End: "; + VLOG(2) << "[3PHS] Bounds: Residual Lower " + << curr_state.dual_state.lower_bound() + model.fixed_cost() + << ", Upper " << best_state.solution.cost(); + + return best_state; +} + +/////////////////////////////////////////////////////////////////////// +///////////////////// OUTER REFINEMENT PROCEDURE ////////////////////// +/////////////////////////////////////////////////////////////////////// + +namespace { + +struct GapContribution { + Cost gap; + FullSubsetIndex idx; +}; + +std::vector SelectColumnByGapContribution( + const SubModel& model, const PrimalDualState& best_state, + BaseInt nrows_to_fix) { + const auto& [solution, dual_state] = best_state; + + FullCoverCounters row_coverage(model.num_elements()); + auto full_model = model.StrongTypedFullModelView(); + + for (FullSubsetIndex j : solution.subsets()) { + row_coverage.Cover(full_model.columns()[j]); + } + + std::vector gap_contributions; + for (FullSubsetIndex j : solution.subsets()) { + Cost j_gap = .0; + Cost reduced_cost = dual_state.reduced_costs()[static_cast(j)]; + for (FullElementIndex i : full_model.columns()[j]) { + Cost i_mult = dual_state.multipliers()[static_cast(i)]; + j_gap += i_mult * (1.0 - 1.0 / row_coverage[i]); + reduced_cost -= i_mult; + } + j_gap += std::max(reduced_cost, 0.0); + gap_contributions.push_back({j_gap, j}); + } + absl::c_sort(gap_contributions, [](GapContribution g1, GapContribution g2) { + return g1.gap < g2.gap; + }); + + BaseInt covered_rows = 0; + row_coverage.Reset(model.num_elements()); + std::vector cols_to_fix; + for (auto [j_gap, j] : gap_contributions) { + covered_rows += row_coverage.Cover(full_model.columns()[j]); + if (covered_rows > nrows_to_fix) { + break; + } + cols_to_fix.push_back(static_cast(j)); + } + return cols_to_fix; +} +} // namespace + +PrimalDualState RunCftHeuristic(SubModel& model, + const Solution& init_solution) { + CFT_MEASURE_SCOPE_DURATION(refinement_time); + + PrimalDualState best_state = {.solution = init_solution, + .dual_state = MakeTentativeDualState(model)}; + if (best_state.solution.Empty()) { + best_state.solution = + RunMultiplierBasedGreedy(model, best_state.dual_state); + } + + Cost cost_cutoff = std::numeric_limits::lowest(); + double fix_minimum = .3; // Arbitrary from [1] + double fix_increment = 1.1; // Arbitrary from [1] + double fix_fraction = fix_minimum; + + for (BaseInt iter_counter = 0; model.num_elements() > 0; ++iter_counter) { + VLOG(1) << "[CFTH] Refinement iteration: " << iter_counter; + Cost fixed_cost_before = model.fixed_cost(); + auto [solution, dual_state] = RunThreePhase(model, best_state.solution); + if (iter_counter == 0) { + best_state.dual_state = std::move(dual_state); + } + if (solution.cost() < best_state.solution.cost()) { + best_state.solution = solution; + fix_fraction = fix_minimum; + } + cost_cutoff = std::max(cost_cutoff, + fixed_cost_before + dual_state.lower_bound()); + VLOG(1) << "[CFTH] Refinement Bounds: Residual Lower " << cost_cutoff + << ", Real Lower " << best_state.dual_state.lower_bound() + << ", Upper " << best_state.solution.cost(); + if (best_state.solution.cost() - CFT_BOUND_EPSILON <= cost_cutoff) { + break; + } + + fix_fraction = std::min(1.0, fix_fraction * fix_increment); + std::vector cols_to_fix = SelectColumnByGapContribution( + model, best_state, model.num_elements() * fix_fraction); + + if (!cols_to_fix.empty()) { + model.ResetColumnFixing(cols_to_fix, best_state.dual_state); + } + VLOG(1) << "[CFTH] Fixed " << cols_to_fix.size() + << " new columns with cost: " << model.fixed_cost(); + VLOG(1) << "[CFTH] Model sizes: rows " << model.num_focus_elements() << "/" + << model.num_elements() << ", columns " << model.num_focus_subsets() + << "/" << model.num_subsets(); + } + +#ifdef CFT_MEASURE_TIME + double subg_t = subgradient_time.Get(); + double greedy_t = greedy_time.Get(); + double three_phase_t = three_phase_time.Get(); + double refinement_t = refinement_time.Get(); + + printf("Subgradient time: %8.2f (%.1f%%)\n", subg_t, + 100 * subg_t / refinement_t); + printf("Greedy Heur time: %8.2f (%.1f%%)\n", greedy_t, + 100 * greedy_t / refinement_t); + printf("SubG - Greedy time: %8.2f (%.1f%%)\n", subg_t - greedy_t, + 100 * (subg_t - greedy_t) / refinement_t); + printf("3Phase time: %8.2f (%.1f%%)\n", three_phase_t, + 100 * three_phase_t / refinement_t); + printf("3Phase - Subg time: %8.2f (%.1f%%)\n", three_phase_t - subg_t, + 100 * (three_phase_t - subg_t) / refinement_t); + printf("Total CFT time: %8.2f (%.1f%%)\n", refinement_t, 100.0); +#endif return best_state; } @@ -760,9 +976,16 @@ PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution) { namespace { std::vector ComputeTentativeFocus(StrongModelView full_model) { - FullSubsetBoolVector selected(full_model.num_subsets(), false); std::vector columns_focus; + + if (full_model.num_subsets() <= 2 * kMinCov * full_model.num_elements()) { + columns_focus.resize(full_model.num_subsets()); + absl::c_iota(columns_focus, FullSubsetIndex(0)); + return columns_focus; + } + columns_focus.reserve(full_model.num_elements() * kMinCov); + FullSubsetBoolVector selected(full_model.num_subsets(), false); // Select the first min_row_coverage columns for each row for (const auto& row : full_model.rows()) { @@ -781,82 +1004,94 @@ std::vector ComputeTentativeFocus(StrongModelView full_model) { absl::c_sort(columns_focus); return columns_focus; } +} // namespace -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); +const std::vector& +FullToCoreModel::ColumnSelector::ComputeNewSelection( + FilterModelView full_model, + const std::vector& forced_columns, + const SubsetCostVector& reduced_costs) { + selected_.assign(full_model.num_subsets(), false); + row_cover_counts_.assign(full_model.num_elements(), 0); + rows_left_to_cover_ = full_model.num_focus_elements(); + selection_.clear(); + candidates_.clear(); + + // Always retain best solution in the core model (if possible) + for (FullSubsetIndex full_j : forced_columns) { + SubsetIndex j = static_cast(full_j); + if (full_model.IsFocusCol(j)) { + SelectColumn(full_model, j); } + } + + auto subset_range = full_model.SubsetRange(); + candidates_ = {subset_range.begin(), subset_range.end()}; + absl::c_sort(candidates_, [&](auto j1, auto j2) { + return reduced_costs[j1] < reduced_costs[j2]; + }); + first_unselected_ = candidates_.begin(); + + SelecteMinRedCostColumns(full_model, reduced_costs); + SelectMinRedCostByRow(full_model, reduced_costs); + absl::c_sort(selection_); + return selection_; +} - 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); +bool FullToCoreModel::ColumnSelector::SelectColumn(FilterModelView full_model, + SubsetIndex j) { + if (selected_[j]) { + return false; + } + for (ElementIndex i : full_model.columns()[j]) { + selected_[j] = true; // Detect empty columns + if (++row_cover_counts_[i] == kMinCov) { + --rows_left_to_cover_; } } + if (selected_[j]) { // Skip empty comlumns + selection_.push_back(static_cast(j)); + } + return selected_[j]; } -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()); +void FullToCoreModel::ColumnSelector::SelecteMinRedCostColumns( + FilterModelView full_model, const SubsetCostVector& reduced_costs) { + BaseInt selected_size = 0; + BaseInt max_size = kMinCov * full_model.num_elements(); + auto it = first_unselected_; + while (it != candidates_.end() && reduced_costs[*it] < 0.1 && + selected_size < max_size) { + selected_size += SelectColumn(full_model, *it++) ? 1 : 0; + } + first_unselected_ = it; +} - 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; - } +void FullToCoreModel::ColumnSelector::SelectMinRedCostByRow( + FilterModelView full_model, const SubsetCostVector& reduced_costs) { + auto it = first_unselected_; + while (it != candidates_.end() && rows_left_to_cover_ > 0) { + SubsetIndex j = *it++; + if (selected_[j]) { + continue; } - - 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); - } + for (ElementIndex i : full_model.columns()[j]) { + ++row_cover_counts_[i]; + rows_left_to_cover_ += (row_cover_counts_[i] == kMinCov ? 1 : 0); + selected_[j] = selected_[j] || (row_cover_counts_[i] <= kMinCov); + } + if (selected_[j]) { + selection_.push_back(static_cast(j)); } } } -} // 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()), + prev_best_lower_bound_(std::numeric_limits::lowest()), full_dual_state_(*full_model), best_dual_state_(full_dual_state_) { ResetPricingPeriod(); @@ -870,15 +1105,19 @@ void FullToCoreModel::ResetPricingPeriod() { update_max_period_ = std::min(1000, full_model_->num_elements() / 3); } -Cost FullToCoreModel::FixColumns( +Cost FullToCoreModel::FixMoreColumns( const std::vector& columns_to_fix) { StrongModelView typed_full_model = StrongTypedFullModelView(); for (SubsetIndex core_j : columns_to_fix) { FullSubsetIndex full_j = SubModel::MapCoreToFullSubsetIndex(core_j); + DCHECK(IsFocusCol(full_j)); IsFocusCol(full_j) = false; + bool any_active = false; for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { + any_active |= IsFocusRow(full_i); IsFocusRow(full_i) = false; } + DCHECK(any_active); } for (FullSubsetIndex full_j : typed_full_model.SubsetRange()) { if (!IsFocusCol(full_j)) { @@ -892,59 +1131,95 @@ Cost FullToCoreModel::FixColumns( } } } + ResetPricingPeriod(); - Cost fixed_cost = base::FixColumns(columns_to_fix); + Cost fixed_cost = base::FixMoreColumns(columns_to_fix); DCHECK(FullToSubModelInvariantCheck()); return fixed_cost; } -bool FullToCoreModel::UpdateCore(PrimalDualState& core_state) { - if (--update_countdown_ > 0) { - return false; +void FullToCoreModel::ResetColumnFixing( + const std::vector& full_columns_to_fix, + const DualState& dual_state) { + is_focus_col_.assign(full_model_->num_subsets(), true); + is_focus_row_.assign(full_model_->num_elements(), true); + + full_dual_state_ = dual_state; + + // We could implement and in-place core-model update that removes old fixings, + // set the new one while also updating the column focus. This solution is much + // simpler. It just create a new core-model object from scratch and then uses + // the existing interface. + const auto& selection = col_selector_.ComputeNewSelection( + FixingFullModelView(), full_columns_to_fix, + full_dual_state_.reduced_costs()); + static_cast(*this) = SubModel(full_model_, selection); + + // TODO(anyone): Improve this. It's Inefficient but hardly a botleneck and it + // also avoid storing a full->core column map. + std::vector columns_to_fix; + for (SubsetIndex core_j : SubsetRange()) { + for (FullSubsetIndex full_j : full_columns_to_fix) { + if (full_j == MapCoreToFullSubsetIndex(core_j)) { + columns_to_fix.push_back(core_j); + break; + } + } } + DCHECK_EQ(columns_to_fix.size(), full_columns_to_fix.size()); + FixMoreColumns(columns_to_fix); + DCHECK(FullToSubModelInvariantCheck()); +} - 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'; +void FullToCoreModel::SizeUpdate() { + is_focus_col_.resize(full_model_->num_subsets(), true); +} - 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; - } +bool FullToCoreModel::IsTimeToUpdate(Cost best_lower_bound, bool force) { + if (!force && --update_countdown_ > 0) { + return false; } + if (best_lower_bound == prev_best_lower_bound_) { + return false; + } + prev_best_lower_bound_ = best_lower_bound; + return true; +} - 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 - }); +void FullToCoreModel::ComputeAndSetFocus(Cost best_lower_bound, + const Solution& best_solution) { + const auto& selection = col_selector_.ComputeNewSelection( + FixingFullModelView(), best_solution.subsets(), + full_dual_state_.reduced_costs()); + base::SetFocus(selection); + UpdatePricingPeriod(full_dual_state_, best_lower_bound, + best_solution.cost() - fixed_cost()); + VLOG(3) << "[F2CU] Core-update: Lower bounds: real " + << full_dual_state_.lower_bound() << ", core " << best_lower_bound + << ", core size: " << num_focus_elements() << "x" + << num_focus_subsets(); +} +bool FullToCoreModel::UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) { + if (!IsTimeToUpdate(best_lower_bound, force)) { + return false; + } + UpdateMultipliers(best_multipliers); + ComputeAndSetFocus(best_lower_bound, best_solution); DCHECK(FullToSubModelInvariantCheck()); return true; } void FullToCoreModel::UpdatePricingPeriod(const DualState& full_dual_state, - const PrimalDualState& core_state) { - DCHECK_GE(core_state.dual_state.lower_bound() + 1e-6, - full_dual_state.lower_bound()); - DCHECK_GE(core_state.solution.cost(), .0); - - Cost delta = - core_state.dual_state.lower_bound() - full_dual_state.lower_bound(); - Cost ratio = DivideIfGE0(delta, core_state.solution.cost()); + Cost core_lower_bound, + Cost core_upper_bound) { + DCHECK_GE(core_lower_bound + 1e-6, full_dual_state.lower_bound()); + DCHECK_GE(core_upper_bound, .0); + + Cost delta = core_lower_bound - full_dual_state.lower_bound(); + Cost ratio = DivideIfGE0(delta, core_upper_bound); if (ratio <= 1e-6) { update_period_ = std::min(update_max_period_, 10 * update_period_); } else if (ratio <= 0.02) { @@ -957,15 +1232,14 @@ void FullToCoreModel::UpdatePricingPeriod(const DualState& full_dual_state, update_countdown_ = update_period_; } -void FullToCoreModel::UpdateDualState(const DualState& core_dual_state, - DualState& full_dual_state, - DualState& best_dual_state) { +Cost FullToCoreModel::UpdateMultipliers( + const ElementCostVector& core_multipliers) { auto fixing_full_model = FixingFullModelView(); full_dual_state_.DualUpdate( fixing_full_model, [&](ElementIndex full_i, Cost& i_mult) { ElementIndex core_i = MapFullToCoreElementIndex(static_cast(full_i)); - i_mult = core_dual_state.multipliers()[core_i]; + i_mult = core_multipliers[core_i]; }); // Here, we simply check if any columns have been fixed, and only update the @@ -983,6 +1257,7 @@ void FullToCoreModel::UpdateDualState(const DualState& core_dual_state, full_dual_state_.lower_bound() > best_dual_state_.lower_bound()) { best_dual_state_ = full_dual_state_; } + return full_dual_state_.lower_bound(); } bool FullToCoreModel::FullToSubModelInvariantCheck() { @@ -1009,8 +1284,50 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { " not found in full model view.\n"); return false; } + + const auto& core_column = sub_model.columns()[core_j]; + if (core_column.begin() == core_column.end()) { + std::cerr << absl::StrCat("Core subset ", core_j, " empty.\n"); + return false; + } + + const auto& full_column = typed_full_model.columns()[full_j]; + if (full_column.begin() == full_column.end()) { + std::cerr << absl::StrCat("Full subset ", full_j, " empty.\n"); + return false; + } + + // Assumes corresponding elements have the same order in both models + auto core_it = core_column.begin(); + for (FullElementIndex full_i : typed_full_model.columns()[full_j]) { + if (sub_model.MapFullToCoreElementIndex(full_i) != *core_it) { + continue; + } + if (sub_model.MapCoreToFullElementIndex(*core_it) != full_i) { + std::cerr << absl::StrCat( + "Subset ", core_j, " in sub-model has mapped element ", *core_it, + " but it is not the same as the full model.\n"); + return false; + } + if (++core_it == core_column.end()) { + break; + } + } + if (core_it != core_column.end()) { + std::cerr << absl::StrCat("Subset ", core_j, + " in sub-model has no mapped element ", + *core_it, " in full model view.\n"); + return false; + } } + for (ElementIndex core_i : sub_model.ElementRange()) { + const auto& core_row = sub_model.rows()[core_i]; + if (core_row.begin() == core_row.end()) { + std::cerr << absl::StrCat("Core row ", core_i, " empty.\n"); + return false; + } + FullElementIndex full_i = sub_model.MapCoreToFullElementIndex(core_i); if (!is_focus_row_[static_cast(full_i)]) { std::cerr << absl::StrCat("Element ", core_i, @@ -1019,6 +1336,7 @@ bool FullToCoreModel::FullToSubModelInvariantCheck() { return false; } } + for (FullElementIndex full_i : typed_full_model.ElementRange()) { if (!IsFocusRow(full_i)) { continue; diff --git a/ortools/set_cover/set_cover_cft.h b/ortools/set_cover/set_cover_cft.h index 5cdc99123eb..9481161a570 100644 --- a/ortools/set_cover/set_cover_cft.h +++ b/ortools/set_cover/set_cover_cft.h @@ -100,7 +100,7 @@ class Solution { } private: - Cost cost_; + Cost cost_ = std::numeric_limits::max(); std::vector subsets_; }; @@ -122,8 +122,9 @@ inline Cost DivideIfGE0(Cost numerator, Cost denominator) { class DualState { public: DualState() = default; + DualState(const DualState&) = default; template - DualState(const SubModelT& model) + explicit DualState(const SubModelT& model) : lower_bound_(.0), multipliers_(model.num_elements(), .0), reduced_costs_(model.subset_costs().begin(), @@ -143,6 +144,7 @@ class DualState { for (ElementIndex i : model.ElementRange()) { multiplier_operator(i, multipliers_[i]); lower_bound_ += multipliers_[i]; + DCHECK(std::isfinite(multipliers_[i])); DCHECK_GE(multipliers_[i], .0); } lower_bound_ += ComputeReducedCosts(model, multipliers_, reduced_costs_); @@ -188,7 +190,11 @@ struct PrimalDualState { struct SubgradientContext { const SubModel& model; const DualState& current_dual_state; - const DualState& best_dual_state; + + // Avoid copying unused reduced cost during subgradient + const Cost& best_lower_bound; + const ElementCostVector& best_multipliers; + const Solution& best_solution; const ElementCostVector& subgradient; }; @@ -201,7 +207,8 @@ class SubgradientCBs { virtual void RunHeuristic(const SubgradientContext&, Solution&) = 0; virtual void ComputeMultipliersDelta(const SubgradientContext&, ElementCostVector& delta_mults) = 0; - virtual bool UpdateCoreModel(SubModel&, PrimalDualState&) = 0; + virtual bool UpdateCoreModel(SubgradientContext context, + CoreModel& core_model, bool force = false) = 0; virtual ~SubgradientCBs() = default; }; @@ -223,8 +230,8 @@ class BoundCBs : public SubgradientCBs { ElementCostVector& delta_mults) override; void RunHeuristic(const SubgradientContext& context, Solution& solution) override {} - bool UpdateCoreModel(SubModel& core_model, - PrimalDualState& best_state) override; + bool UpdateCoreModel(SubgradientContext context, CoreModel& core_model, + bool force = false) override; private: void MakeMinimalCoverageSubgradient(const SubgradientContext& context, @@ -232,24 +239,7 @@ class BoundCBs : public SubgradientCBs { private: Cost squared_norm_; - - // This addition implements a simplified stabilization technique inspired by - // works in the literature, such as: - // Frangioni, A., Gendron, B., Gorgone, E. (2015): - // "On the computational efficiency of subgradient methods: A case study in - // combinatorial optimization." - // - // The approach aims to reduce oscillations (zig-zagging) in the subgradient - // by using a "moving average" of the current and previous subgradients. The - // current subgradient is weighted by a factor of alpha, while the previous - // subgradients contribution is weighted by (1 - alpha). The parameter alpha - // is set to 0.5 by default but can be adjusted for tuning. The resulting - // stabilized subgradient is referred to as "direction" to distinguish it from - // the original subgradient. - Cost stabilization_coeff = 0.5; // Arbitrary from c4v4 ElementCostVector direction_; - ElementCostVector prev_direction_; - std::vector lagrangian_solution_; // Stopping condition @@ -257,10 +247,10 @@ class BoundCBs : public SubgradientCBs { BaseInt max_iter_countdown_; BaseInt exit_test_countdown_; BaseInt exit_test_period_; - BaseInt last_core_update_countdown_; + BaseInt unfixed_run_extension_; // Step size - void UpdateStepSize(Cost lower_bound); + void UpdateStepSize(SubgradientContext context); Cost step_size_; Cost last_min_lb_seen_; Cost last_max_lb_seen_; @@ -301,14 +291,15 @@ class HeuristicCBs : public SubgradientCBs { bool ExitCondition(const SubgradientContext& context) override { Cost upper_bound = context.best_solution.cost() - context.model.fixed_cost(); - Cost lower_bound = context.best_dual_state.lower_bound(); + Cost lower_bound = context.best_lower_bound; return upper_bound - .999 < lower_bound || --countdown_ <= 0; } void RunHeuristic(const SubgradientContext& context, Solution& solution) override; void ComputeMultipliersDelta(const SubgradientContext& context, ElementCostVector& delta_mults) override; - bool UpdateCoreModel(SubModel& model, PrimalDualState& state) override { + bool UpdateCoreModel(SubgradientContext context, CoreModel& core_model, + bool force = false) override { return false; } @@ -320,10 +311,20 @@ class HeuristicCBs : public SubgradientCBs { PrimalDualState RunThreePhase(SubModel& model, const Solution& init_solution = {}); +/////////////////////////////////////////////////////////////////////// +///////////////////// OUTER REFINEMENT PROCEDURE ////////////////////// +/////////////////////////////////////////////////////////////////////// + +PrimalDualState RunCftHeuristic(SubModel& model, + const Solution& init_solution = {}); + /////////////////////////////////////////////////////////////////////// //////////////////////// FULL TO CORE PRICING ///////////////////////// /////////////////////////////////////////////////////////////////////// +// Coverage counter to decide the number of columns to keep in the core model. +static constexpr BaseInt kMinCov = 5; + // CoreModel extractor. Stores a pointer to the full model and specilized // UpdateCore in such a way to updated the SubModel (stored as base class) and // focus the search on a small windows of the full model. @@ -335,14 +336,63 @@ class FullToCoreModel : public SubModel { BaseInt max_period; }; + // This class handles the logic for selecting columns based on their reduced + // costs and the number of rows they cover. While this implementation is more + // complex than what would typically be required for static `SetCoverModel`s, + // it is designed to efficiently handle dynamically updated models where new + // columns are generated over time. In this case, recomputing the row view + // from scratch each time would introduce significant overhead. To avoid this, + // the column selection logic operates solely on the column view, without + // relying on the row view. + // + // NOTE: A cleaner alternative would involve modifying the `SetCoverModel` + // implementation to support incremental updates to the row view as new + // columns are added. This approach would reduce overhead while enabling a + // simpler and more efficient column selection process. + // + // NOTE: The row-view based approach is available at commit: + // a598cf83930629853f72b964ebcff01f7a9378e0 + class ColumnSelector { + public: + const std::vector& ComputeNewSelection( + FilterModelView full_model, + const std::vector& forced_columns, + const SubsetCostVector& reduced_costs); + + private: + bool SelectColumn(FilterModelView full_model, SubsetIndex j); + void SelecteMinRedCostColumns(FilterModelView full_model, + const SubsetCostVector& reduced_costs); + void SelectMinRedCostByRow(FilterModelView full_model, + const SubsetCostVector& reduced_costs); + + private: + std::vector candidates_; + std::vector::const_iterator first_unselected_; + ElementToIntVector row_cover_counts_; + BaseInt rows_left_to_cover_; + + std::vector selection_; + SubsetBoolVector selected_; + }; + public: + FullToCoreModel() = default; FullToCoreModel(const Model* full_model); - Cost FixColumns(const std::vector& columns_to_fix) override; - bool UpdateCore(PrimalDualState& core_state) override; + Cost FixMoreColumns(const std::vector& columns_to_fix) override; + void ResetColumnFixing(const std::vector& columns_to_fix, + const DualState& state) override; + bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) override; void ResetPricingPeriod(); const DualState& best_dual_state() const { return best_dual_state_; } + bool FullToSubModelInvariantCheck(); + + protected: + void SizeUpdate(); + bool IsTimeToUpdate(Cost best_lower_bound, bool force); - private: decltype(auto) IsFocusCol(FullSubsetIndex j) { return is_focus_col_[static_cast(j)]; } @@ -350,9 +400,9 @@ class FullToCoreModel : public SubModel { return is_focus_row_[static_cast(i)]; } void UpdatePricingPeriod(const DualState& full_dual_state, - const PrimalDualState& core_state); - void UpdateDualState(const DualState& core_dual_state, - DualState& full_dual_state, DualState& best_dual_state); + Cost core_lower_bound, Cost core_upper_bound); + Cost UpdateMultipliers(const ElementCostVector& core_multipliers); + void ComputeAndSetFocus(Cost best_lower_bound, const Solution& best_solution); // Views are not composable (for now), so we can either access the full_model // with the strongly typed view or with the filtered view. @@ -360,7 +410,8 @@ class FullToCoreModel : public SubModel { // Access the full model filtered by the current columns fixed. FilterModelView FixingFullModelView() const { return FilterModelView(full_model_, &is_focus_col_, &is_focus_row_, - num_subsets_, num_elements_); + full_model_->num_subsets(), + full_model_->num_elements()); } // Access the full model with the strongly typed view. @@ -368,7 +419,8 @@ class FullToCoreModel : public SubModel { return StrongModelView(full_model_); } - bool FullToSubModelInvariantCheck(); + std::vector SelectNewCoreColumns( + const std::vector& forced_columns = {}); private: const Model* full_model_; @@ -388,19 +440,17 @@ class FullToCoreModel : public SubModel { // implementation that avoids this memory optimization was preferred. ElementBoolVector is_focus_row_; - BaseInt num_subsets_; - BaseInt num_elements_; - + BaseInt selection_coefficient_ = kMinCov; + Cost prev_best_lower_bound_; DualState full_dual_state_; DualState best_dual_state_; BaseInt update_countdown_; BaseInt update_period_; BaseInt update_max_period_; -}; -// Coverage counter to decide the number of columns to keep in the core model. -static constexpr BaseInt kMinCov = 5; + ColumnSelector col_selector_; // Here to avoid reallocations +}; } // namespace operations_research::scp diff --git a/ortools/set_cover/set_cover_submodel.cc b/ortools/set_cover/set_cover_submodel.cc index 02a13eb2559..59202e9f925 100644 --- a/ortools/set_cover/set_cover_submodel.cc +++ b/ortools/set_cover/set_cover_submodel.cc @@ -14,37 +14,15 @@ #include "ortools/set_cover/set_cover_submodel.h" #include "ortools/base/stl_util.h" +#include "ortools/set_cover/set_cover_cft.h" +#include "ortools/set_cover/set_cover_views.h" namespace operations_research::scp { -template -void PrintSubModelSummary(const SubModelT& model) { - std::cout << "SubModelView sizes: rows " << model.num_focus_elements() << "/" - << model.num_elements() << ", columns " << model.num_focus_subsets() - << "/" << model.num_subsets() << '\n'; - std::cout << "Fixing: columns " << model.fixed_columns().size() << " cost " - << model.fixed_cost() << '\n'; -} - SubModelView::SubModelView(const Model* model) : base_view(model, &cols_sizes_, &rows_sizes_, &cols_focus_, &rows_focus_), - full_model_(model), - cols_sizes_(model->num_subsets()), - rows_sizes_(model->num_elements()) { - DCHECK(full_model_ != nullptr); - cols_focus_.reserve(model->num_subsets()); - rows_focus_.reserve(model->num_elements()); - for (SubsetIndex j : model->SubsetRange()) { - cols_sizes_[j] = model->columns()[j].size(); - cols_focus_.push_back(j); - } - for (ElementIndex i : model->ElementRange()) { - rows_sizes_[i] = model->rows()[i].size(); - rows_focus_.push_back(i); - } - fixed_columns_.clear(); - fixed_cost_ = 0.0; - PrintSubModelSummary(*this); + full_model_(model) { + ResetToIdentitySubModel(); DCHECK(ValidateSubModel(*this)); } @@ -59,7 +37,25 @@ SubModelView::SubModelView(const Model* model, SetFocus(columns_focus); } -Cost SubModelView::FixColumns(const std::vector& columns_to_fix) { +void SubModelView::ResetToIdentitySubModel() { + cols_sizes_.resize(full_model_->num_subsets()); + rows_sizes_.resize(full_model_->num_elements()); + cols_focus_.clear(); + rows_focus_.clear(); + for (SubsetIndex j : full_model_->SubsetRange()) { + cols_sizes_[j] = full_model_->columns()[j].size(); + cols_focus_.push_back(j); + } + for (ElementIndex i : full_model_->ElementRange()) { + rows_sizes_[i] = full_model_->rows()[i].size(); + rows_focus_.push_back(i); + } + fixed_columns_.clear(); + fixed_cost_ = .0; +} + +Cost SubModelView::FixMoreColumns( + const std::vector& columns_to_fix) { DCHECK(full_model_ != nullptr); Cost old_fixed_cost = fixed_cost_; if (columns_to_fix.empty()) { @@ -88,11 +84,20 @@ Cost SubModelView::FixColumns(const std::vector& columns_to_fix) { gtl::STLEraseAllFromSequenceIf( &rows_focus_, [&](ElementIndex i) { return !rows_sizes_[i]; }); - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); return fixed_cost_ - old_fixed_cost; } +void SubModelView::ResetColumnFixing( + const std::vector& columns_to_fix, const DualState&) { + ResetToIdentitySubModel(); + std::vector core_column_to_fix; + for (FullSubsetIndex full_j : columns_to_fix) { + core_column_to_fix.push_back(static_cast(full_j)); + } + FixMoreColumns(core_column_to_fix); +} + void SubModelView::SetFocus(const std::vector& columns_focus) { DCHECK(full_model_ != nullptr); DCHECK(!rows_sizes_.empty()); @@ -125,24 +130,15 @@ void SubModelView::SetFocus(const std::vector& columns_focus) { rows_focus_.push_back(i); } } - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); } -CoreModel::CoreModel(const Model* model) - : Model(*model), - full_model_(model), - core2full_row_map_(model->num_elements()), - full2core_row_map_(model->num_elements()), - core2full_col_map_(model->num_subsets()) { +CoreModel::CoreModel(const Model* model) : Model(), full_model_(model) { CHECK(ElementIndex(full_model_.num_elements()) < null_element_index) << "Max element index is reserved."; CHECK(SubsetIndex(full_model_.num_subsets()) < null_subset_index) << "Max subset index is reserved."; - - absl::c_iota(core2full_row_map_, FullElementIndex()); - absl::c_iota(full2core_row_map_, ElementIndex()); - absl::c_iota(core2full_col_map_, FullSubsetIndex()); + ResetToIdentitySubModel(); } CoreModel::CoreModel(const Model* model, @@ -161,6 +157,18 @@ CoreModel::CoreModel(const Model* model, SetFocus(columns_focus); } +void CoreModel::ResetToIdentitySubModel() { + core2full_row_map_.resize(full_model_.num_elements()); + full2core_row_map_.resize(full_model_.num_elements()); + core2full_col_map_.resize(full_model_.num_subsets()); + absl::c_iota(core2full_row_map_, FullElementIndex()); + absl::c_iota(full2core_row_map_, ElementIndex()); + absl::c_iota(core2full_col_map_, FullSubsetIndex()); + fixed_cost_ = .0; + fixed_columns_.clear(); + static_cast(*this) = Model(full_model_.base()); +} + // Note: Assumes that columns_focus covers all rows for which rows_flags is true // (i.e.: non-covered rows should be set to false to rows_flags). This property // get exploited to keep the rows in the same ordering of the original model @@ -177,7 +185,6 @@ void CoreModel::SetFocus(const std::vector& columns_focus) { // Now we can fill the new core model for (FullSubsetIndex full_j : columns_focus) { - core2full_col_map_.push_back(full_j); bool first_row = true; for (FullElementIndex full_i : full_model_.columns()[full_j]) { ElementIndex core_i = full2core_row_map_[full_i]; @@ -190,10 +197,13 @@ void CoreModel::SetFocus(const std::vector& columns_focus) { submodel.AddElementToLastSubset(core_i); } } + // Handle empty columns + if (!first_row) { + core2full_col_map_.push_back(full_j); + } } submodel.CreateSparseRowView(); - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); } @@ -273,7 +283,7 @@ Model CoreModel::MakeNewCoreModel( return new_submodel; } -Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { +Cost CoreModel::FixMoreColumns(const std::vector& columns_to_fix) { if (columns_to_fix.empty()) { return .0; } @@ -288,7 +298,6 @@ Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { // Create new model object applying the computed mappings static_cast(*this) = MakeNewCoreModel(new_c2f_row_map); - PrintSubModelSummary(*this); DCHECK(ValidateSubModel(*this)); DCHECK(absl::c_is_sorted(core2full_col_map_)); DCHECK(absl::c_is_sorted(core2full_row_map_)); @@ -296,4 +305,14 @@ Cost CoreModel::FixColumns(const std::vector& columns_to_fix) { return fixed_cost_ - old_fixed_cost; } +void CoreModel::ResetColumnFixing( + const std::vector& columns_to_fix, const DualState&) { + ResetToIdentitySubModel(); + std::vector core_column_to_fix; + for (FullSubsetIndex full_j : columns_to_fix) { + core_column_to_fix.push_back(static_cast(full_j)); + } + FixMoreColumns(core_column_to_fix); +} + } // namespace operations_research::scp \ No newline at end of file diff --git a/ortools/set_cover/set_cover_submodel.h b/ortools/set_cover/set_cover_submodel.h index 18dc20d29e9..13b4ce93f1c 100644 --- a/ortools/set_cover/set_cover_submodel.h +++ b/ortools/set_cover/set_cover_submodel.h @@ -26,6 +26,7 @@ using Model = SetCoverModel; // Forward declarations, see below for the definition of the classes. struct PrimalDualState; struct Solution; +struct DualState; // The CFT algorithm generates sub-models in two distinct ways: // @@ -112,15 +113,29 @@ class SubModelView : public IndexListModelView { // Fix the provided columns, removing them for the submodel. Rows now covered // by fixed columns are also removed from the submodel along with non-fixed // columns that only cover those rows. - virtual Cost FixColumns(const std::vector& columns_to_fix); + virtual Cost FixMoreColumns(const std::vector& columns_to_fix); + + virtual void ResetColumnFixing( + const std::vector& columns_to_fix, + const DualState& state); // Hook function for specializations. This function can be used to define a // "small" core model considering a subset of the full model through the use // of column-generation or by only selecting columns with good reduced cost in // the full model. - virtual bool UpdateCore(PrimalDualState& core_state) { return false; } + virtual bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) { + return false; + } + + StrongModelView StrongTypedFullModelView() const { + return StrongModelView(full_model_); + } private: + void ResetToIdentitySubModel(); + // Pointer to the original model const Model* full_model_; @@ -179,7 +194,6 @@ class CoreModel : private Model { ElementIndex MapFullToCoreElementIndex(FullElementIndex full_i) const { DCHECK(FullElementIndex() <= full_i && full_i < FullElementIndex(num_elements())); - DCHECK(full2core_row_map_[full_i] != null_element_index); return full2core_row_map_[full_i]; } FullSubsetIndex MapCoreToFullSubsetIndex(SubsetIndex core_j) const { @@ -210,18 +224,29 @@ class CoreModel : private Model { // Fix the provided columns, removing them for the submodel. Rows now covered // by fixed columns are also removed from the submodel along with non-fixed // columns that only cover those rows. - virtual Cost FixColumns(const std::vector& columns_to_fix); + virtual Cost FixMoreColumns(const std::vector& columns_to_fix); + + virtual void ResetColumnFixing( + const std::vector& columns_to_fix, + const DualState& state); // Hook function for specializations. This function can be used to define a // "small" core model considering a subset of the full model through the use // of column-generation or by only selecting columns with good reduced cost in // the full model. - virtual bool UpdateCore(PrimalDualState& core_state) { return false; } + virtual bool UpdateCore(Cost best_lower_bound, + const ElementCostVector& best_multipliers, + const Solution& best_solution, bool force) { + return false; + } + + StrongModelView StrongTypedFullModelView() const { return full_model_; } private: void MarkNewFixingInMaps(const std::vector& columns_to_fix); CoreToFullElementMapVector MakeOrFillBothRowMaps(); Model MakeNewCoreModel(const CoreToFullElementMapVector& new_c2f_col_map); + void ResetToIdentitySubModel(); // Pointer to the original model StrongModelView full_model_; diff --git a/ortools/set_cover/set_cover_views.h b/ortools/set_cover/set_cover_views.h index d9c97d725c1..12e30390ebc 100644 --- a/ortools/set_cover/set_cover_views.h +++ b/ortools/set_cover/set_cover_views.h @@ -57,6 +57,10 @@ using FullElementCostVector = util_intops::StrongVector; using FullSubsetCostVector = util_intops::StrongVector; using FullElementBoolVector = util_intops::StrongVector; using FullSubsetBoolVector = util_intops::StrongVector; +using FullElementToIntVector = + util_intops::StrongVector; +using FullSubsetToIntVector = + util_intops::StrongVector; // When a sub-model is created, indicies are compacted to be consecutive and // strarting from 0 (to reduce memory usage). Core ElementIndex to original @@ -121,6 +125,7 @@ class StrongModelView { auto ElementRange() const -> util_intops::StrongIntRange { return {FullElementIndex(), FullElementIndex(num_elements())}; } + const SetCoverModel& base() const { return *model_; } private: const SetCoverModel* model_; @@ -180,6 +185,7 @@ class IndexListModelView { DCHECK(ElementIndex() <= i && i < ElementIndex(num_elements())); return (*rows_sizes_)[i]; } + const SetCoverModel& base() const { return *model_; } private: const SetCoverModel* model_; @@ -233,6 +239,10 @@ class FilterModelView { -> util_intops::FilterIndexRangeView { return {is_focus_row_}; } + bool IsFocusCol(SubsetIndex j) const { return (*is_focus_col_)[j]; } + bool IsFocusRow(ElementIndex i) const { return (*is_focus_row_)[i]; } + + const SetCoverModel& base() const { return *model_; } private: const SetCoverModel* model_; diff --git a/ortools/set_cover/views.h b/ortools/set_cover/views.h index ba4c5a46c38..49a31be0b02 100644 --- a/ortools/set_cover/views.h +++ b/ortools/set_cover/views.h @@ -19,6 +19,8 @@ #include +#include "ortools/graph/iterators.h" + // NOTE: It may be unexpected, but views provide a subscript operator that // directly accesses the underlying original container using the original // indices. This behavior is particularly relevant in the context of the Set @@ -180,6 +182,7 @@ class IndexListView { IndexListViewIterator end() const { return IndexListViewIterator(values_, indices_.end()); } + absl::Span base() const { return values_; } private: absl::Span values_; @@ -255,6 +258,7 @@ class ValueFilterView { << "Inactive value. Are you using relative indices?"; return value; } + absl::Span base() const { return values_; } private: absl::Span values_; @@ -329,6 +333,7 @@ class IndexFilterView { return IndexFilterViewIterator(values_.end(), is_active_->end(), is_active_->end()); } + absl::Span base() const { return values_; } // NOTE: uses indices of the original container, not the filtered one template @@ -374,6 +379,7 @@ class TwoLevelsView : public Lvl1ViewT { level2_type operator*() const { return level2_type(&(*iter_), active_items_); } + const Lvl1ViewT& base() const { return *this; } private: level1_iterator iter_; @@ -460,6 +466,7 @@ class TransformView { TransformViewIterator end() const { return TransformViewIterator(values_.end()); } + absl::Span base() const { return values_; } private: absl::Span values_;