22#include "absl/log/check.h"
23#include "absl/synchronization/blocking_counter.h"
50 std::numeric_limits<Cost>::infinity());
54 marginal_costs[subset] =
60 Cost min_marginal_cost = std::numeric_limits<Cost>::infinity();
64 for (
const SubsetIndex subset : rows[element]) {
65 min_marginal_cost = std::min(min_marginal_cost, marginal_costs[subset]);
67 multipliers[element] = min_marginal_cost;
79 for (
const ColumnEntryIndex pos : column.
index_range()) {
80 result += dual[column[pos]];
91void FillReducedCostsSlice(SubsetIndex slice_start, SubsetIndex slice_end,
96 for (SubsetIndex subset = slice_start; subset < slice_end; ++subset) {
97 (*reduced_costs)[subset] =
104 return (size + num_threads - 1) / num_threads;
116 const SubsetIndex block_size(BlockSize(
num_subsets.value(), num_threads_));
117 absl::BlockingCounter num_threads_running(num_threads_);
118 SubsetIndex slice_start(0);
119 for (
int thread_index = 0; thread_index < num_threads_; ++thread_index) {
120 const SubsetIndex slice_end =
122 thread_pool_->Schedule([&num_threads_running, slice_start, slice_end,
123 &costs, &multipliers, &columns, &reduced_costs]() {
124 FillReducedCostsSlice(slice_start, slice_end, costs, multipliers, columns,
126 num_threads_running.DecrementCount();
128 slice_start = slice_end;
130 num_threads_running.Wait();
131 return reduced_costs;
143 FillReducedCostsSlice(SubsetIndex(0), SubsetIndex(reduced_costs.size()),
144 costs, multipliers, columns, &reduced_costs);
145 return reduced_costs;
153void FillSubgradientSlice(SubsetIndex slice_start, SubsetIndex slice_end,
157 for (SubsetIndex subset(slice_start); subset < slice_end; ++subset) {
158 if (reduced_costs[subset] < 0.0) {
159 for (
const ElementIndex element : columns[subset]) {
160 (*subgradient)[element] -= 1.0;
176 FillSubgradientSlice(SubsetIndex(0), SubsetIndex(reduced_costs.size()),
177 model()->columns(), reduced_costs, &subgradient);
190 std::vector<ElementCostVector> subgradients(
192 absl::BlockingCounter num_threads_running(num_threads_);
193 const SubsetIndex block_size(BlockSize(
num_subsets.value(), num_threads_));
194 SubsetIndex slice_start(0);
195 for (
int thread_index = 0; thread_index < num_threads_; ++thread_index) {
196 const SubsetIndex slice_end =
198 thread_pool_->Schedule([&num_threads_running, slice_start, slice_end,
199 &reduced_costs, &columns, &subgradients,
201 FillSubgradientSlice(slice_start, slice_end, columns, reduced_costs,
202 &subgradients[thread_index]);
203 num_threads_running.DecrementCount();
205 slice_start = slice_end;
207 num_threads_running.Wait();
208 for (
int thread_index = 0; thread_index < num_threads_; ++thread_index) {
210 subgradient[element] += subgradients[thread_index][element];
223void FillLagrangianValueSlice(SubsetIndex slice_start, SubsetIndex slice_end,
225 Cost* lagrangian_value) {
229 for (SubsetIndex subset(slice_start); subset < slice_end; ++subset) {
230 if (reduced_costs[subset] < 0.0) {
231 *lagrangian_value += reduced_costs[subset];
248 Cost lagrangian_value = 0.0;
250 for (
const Cost u : multipliers) {
251 lagrangian_value += u;
253 FillLagrangianValueSlice(SubsetIndex(0), SubsetIndex(reduced_costs.size()),
254 reduced_costs, &lagrangian_value);
255 return lagrangian_value;
261 Cost lagrangian_value = 0.0;
263 for (
const Cost u : multipliers) {
264 lagrangian_value += u;
266 std::vector<Cost> lagrangian_values(num_threads_, 0.0);
267 absl::BlockingCounter num_threads_running(num_threads_);
268 const SubsetIndex block_size(BlockSize(
model()->
num_subsets(), num_threads_));
270 SubsetIndex slice_start(0);
271 for (
int thread_index = 0; thread_index < num_threads_; ++thread_index) {
272 const SubsetIndex slice_end =
274 thread_pool_->Schedule([&num_threads_running, slice_start, block_size,
276 &lagrangian_values]() {
277 const SubsetIndex slice_end =
279 FillLagrangianValueSlice(slice_start, slice_end, reduced_costs,
280 &lagrangian_values[thread_index]);
281 num_threads_running.DecrementCount();
283 slice_start = slice_end;
285 num_threads_running.Wait();
286 for (
const Cost l : lagrangian_values) {
287 lagrangian_value += l;
289 return lagrangian_value;
307 double step_size,
Cost lagrangian_value,
Cost upper_bound,
311 DCHECK_GT(step_size, 0);
314 Cost subgradient_square_norm = 0.0;
315 for (
const Cost x : subgradient) {
316 subgradient_square_norm += x * x;
320 step_size * (upper_bound - lagrangian_value) / subgradient_square_norm;
324 (*multipliers)[element] = std::clamp(
325 (*multipliers)[element] + factor * subgradient[element], 0.0, 1e6);
330 double step_size,
Cost lagrangian_value,
Cost upper_bound,
334 DCHECK_GT(step_size, 0);
338 Cost subgradient_square_norm = 0.0;
339 for (
const Cost x : subgradient) {
340 subgradient_square_norm += x * x;
344 step_size * (upper_bound - lagrangian_value) / subgradient_square_norm;
347 const Cost kRoof = 1e6;
348 (*multipliers)[element] = std::clamp(
349 (*multipliers)[element] + factor * subgradient[element], 0.0, kRoof);
359 if (
solution[subset] && reduced_costs[subset] > 0.0) {
360 gap += reduced_costs[subset];
361 }
else if (!
solution[subset] && reduced_costs[subset] < 0.0) {
363 gap -= reduced_costs[subset];
368 gap += (coverage[element] - 1) * multipliers[element];
382 delta[subset] = std::max(reduced_costs[subset], 0.0);
383 for (
const ElementIndex element : columns[subset]) {
384 const double size = coverage[element];
385 delta[subset] += multipliers[element] * (size - 1.0) / size;
400 StepSizer(
int period,
double step_size)
401 : period_(period), iter_to_check_(period), step_size_(step_size) {
404 double UpdateStepSize(
int iter,
Cost lower_bound) {
405 min_lb_ = std::min(min_lb_, lower_bound);
406 max_lb_ = std::max(max_lb_, lower_bound);
407 if (iter == iter_to_check_) {
408 iter_to_check_ += period_;
413 max_lb_ == 0.0 ? 0.0 : (max_lb_ - min_lb_) / std::abs(max_lb_);
414 DCHECK_GE(lb_gap, 0.0);
417 }
else if (lb_gap <= 0.001) {
420 step_size_ = std::clamp(step_size_, 1e-6, 10.0);
428 min_lb_ = std::numeric_limits<Cost>::infinity();
429 max_lb_ = -std::numeric_limits<Cost>::infinity();
443 explicit Stopper(
int period)
445 iter_to_check_(period),
446 lower_bound_(std::numeric_limits<
Cost>::min()) {}
447 bool DecideWhetherToStop(
int iter,
Cost lb) {
448 DCHECK_GE(lb, lower_bound_);
449 if (iter == iter_to_check_) {
450 iter_to_check_ += period_;
451 const Cost delta = lb - lower_bound_;
452 const Cost relative_delta = delta / lb;
453 DCHECK_GE(delta, 0.0);
454 DCHECK_GE(relative_delta, 0.0);
456 return relative_delta < 0.001 && delta < 1;
471template <
typename Priority,
typename Index,
bool IsMaxHeap>
474 explicit TopKHeap(
int max_size) : heap_(), max_size_(max_size) {}
475 void Clear() { heap_.Clear(); }
476 void Add(Priority priority, Index index) {
477 if (heap_.Size() < max_size_) {
478 heap_.Add(priority, index);
479 }
else if (heap_.HasPriority(priority, heap_.TopPriority())) {
480 heap_.ReplaceTop(priority, index);
485 AdjustableKAryHeap<Priority,
Index, 2, !IsMaxHeap> heap_;
491std::tuple<Cost, SubsetCostVector, ElementCostVector>
495 Cost lower_bound = 0.0;
497 double step_size = 0.1;
498 StepSizer step_sizer(20, step_size);
499 Stopper stopper(100);
505 for (
int iter = 0; iter < 1000; ++iter) {
507 const Cost lagrangian_value =
510 reduced_costs, &multipliers);
511 lower_bound = std::max(lower_bound, lagrangian_value);
520 return std::make_tuple(lower_bound, reduced_costs, multipliers);
void ReportLowerBound(Cost lower_bound, bool is_cost_consistent)
const ElementToIntVector & coverage() const
Returns vector containing number of subsets covering each element.
void UpdateMultipliers(double step_size, Cost lagrangian_value, Cost upper_bound, const SubsetCostVector &reduced_costs, ElementCostVector *multipliers) const
Cost ComputeLagrangianValue(const SubsetCostVector &reduced_costs, const ElementCostVector &multipliers) const
ElementCostVector ComputeSubgradient(const SubsetCostVector &reduced_costs) const
void ParallelUpdateMultipliers(double step_size, Cost lagrangian_value, Cost upper_bound, const SubsetCostVector &reduced_costs, ElementCostVector *multipliers) const
ElementCostVector ParallelComputeSubgradient(const SubsetCostVector &reduced_costs) const
SubsetCostVector ParallelComputeReducedCosts(const SubsetCostVector &costs, const ElementCostVector &multipliers) const
Computes the reduced costs for all subsets in parallel using ThreadPool.
Cost ParallelComputeLagrangianValue(const SubsetCostVector &reduced_costs, const ElementCostVector &multipliers) const
ElementCostVector InitializeLagrangeMultipliers() const
Initializes the multipliers vector (u) based on the cost per subset.
std::tuple< Cost, SubsetCostVector, ElementCostVector > ComputeLowerBound(const SubsetCostVector &costs, Cost upper_bound)
Computes a lower bound on the optimal cost.
SubsetCostVector ComputeReducedCosts(const SubsetCostVector &costs, const ElementCostVector &multipliers) const
Cost ComputeGap(const SubsetCostVector &reduced_costs, const SubsetBoolVector &solution, const ElementCostVector &multipliers) const
const SparseRowView & rows() const
Row view of the set covering problem.
const SubsetCostVector & subset_costs() const
Vector of costs for each subset.
const SparseColumnView & columns() const
Column view of the set covering problem.
BaseInt num_subsets() const
absl::Duration run_time_
run_time_ is an abstract duration for the time spent in NextSolution().
SetCoverModel * model() const
Accessors.
SetCoverInvariant * inv() const
StrongIntRange< IntType > index_range() const
Fractional ScalarProduct(const DenseRowOrColumn1 &u, const DenseRowOrColumn2 &v)
In SWIG mode, we don't want anything besides these top-level includes.
util_intops::StrongVector< SubsetIndex, Cost > SubsetCostVector
Select next search node to expand Select next item_i to add this new search node to the search Generate a new search node where item_i is not in the knapsack Check validity of this new partial solution(using propagators) - If valid
util_intops::StrongIntRange< SubsetIndex > SubsetRange
util_intops::StrongVector< ElementIndex, Cost > ElementCostVector
util_intops::StrongVector< ElementIndex, BaseInt > ElementToIntVector
util_intops::StrongVector< SubsetIndex, bool > SubsetBoolVector
util_intops::StrongVector< SubsetIndex, SparseColumn > SparseColumnView
Views of the sparse vectors.
double Cost
Basic non-strict type for cost. The speed penalty for using double is ~2%.
util_intops::StrongIntRange< ElementIndex > ElementRange
util_intops::StrongVector< ColumnEntryIndex, ElementIndex > SparseColumn
util_intops::StrongVector< ElementIndex, SparseRow > SparseRowView