set_cover: backport from main

This commit is contained in:
Corentin Le Molgat
2025-04-30 15:12:44 +02:00
parent 8480ff4f3f
commit 6351d3b3f1
21 changed files with 3640 additions and 200 deletions

View File

@@ -2,69 +2,6 @@
This directory contains data structures and algorithms for various problems.
## Set covering
An instance of set covering is composed of two entities: elements and sets; sets
cover a series of elements. The problem of set covering is about finding the
cheapest combination of sets that cover all the elements.
[More information on Wikipedia](https://en.wikipedia.org/wiki/Set_cover_problem).
* Solver: [`set_cover_heuristics.h`](set_cover_heuristics.h).
* Instance representation: [`set_cover_model.h`](set_cover_model.h).
* Instance parser: [`set_cover_reader.h`](set_cover_reader.h).
Create an instance:
```cpp
// If the elements are integers, a subset can be a std::vector<int> (in a pair
// along its cost).
std::vector<std::pair<std::vector<int>, int>> subsets = ...;
SetCoverModel model;
for (const auto [subset, subset_cost] : subsets) {
model.AddEmptySubset(subset_cost)
for (const int element : subset) {
model.AddElementToLastSubset(element);
}
}
SetCoverLedger ledger(&model);
```
Solve it using a MIP solver (guarantees to yield the optimum solution if it has
enough time to run):
```cpp
SetCoverMip mip(&ledger);
mip.SetTimeLimitInSeconds(10);
mip.NextSolution();
SubsetBoolVector best_choices = ledger.GetSolution();
LOG(INFO) << "Cost: " << ledger.cost();
```
A custom combination of heuristics (10,000 iterations of: clearing 10% of the
variables, running a Chvatal greedy descent, using steepest local search):
```cpp
Cost best_cost = std::numeric_limits<Cost>::max();
SubsetBoolVector best_choices = ledger.GetSolution();
for (int i = 0; i < 10000; ++i) {
ledger.LoadSolution(best_choices);
ClearRandomSubsets(0.1 * model.num_subsets().value(), &ledger);
GreedySolutionGenerator greedy(&ledger);
CHECK(greedy.NextSolution());
SteepestSearch steepest(&ledger);
CHECK(steepest.NextSolution(10000));
EXPECT_TRUE(ledger.CheckSolution());
if (ledger.cost() < best_cost) {
best_cost = ledger.cost();
best_choices = ledger.GetSolution();
LOG(INFO) << "Better cost: " << best_cost << " at iteration = " << i;
}
}
ledger.LoadSolution(best_choices);
LOG(INFO) << "Best cost: " << ledger.cost();
```
It used to contain a solver for set covering, which has moved to the
[set_cover](../set_cover)
folder.

View File

@@ -46,7 +46,20 @@ cc_library(
"//ortools/base:intops",
"//ortools/base:strong_vector",
"//ortools/base:timer",
"@abseil-cpp//absl/log:check",
"@abseil-cpp//absl/time",
"@abseil-cpp//absl/types:span",
],
)
cc_test(
name = "base_types_test",
srcs = ["base_types_test.cc"],
deps = [
":base_types",
"//ortools/base:gmock_main",
"@abseil-cpp//absl/types:span",
"@com_google_benchmark//:benchmark",
],
)
@@ -62,6 +75,7 @@ cc_library(
"//ortools/algorithms:adjustable_k_ary_heap",
"//ortools/base:threadpool",
"@abseil-cpp//absl/log:check",
"@abseil-cpp//absl/strings:string_view",
"@abseil-cpp//absl/synchronization",
"@abseil-cpp//absl/types:span",
],
@@ -84,6 +98,7 @@ cc_library(
"@abseil-cpp//absl/random:distributions",
"@abseil-cpp//absl/strings",
"@abseil-cpp//absl/strings:str_format",
"@abseil-cpp//absl/time",
"@abseil-cpp//absl/types:span",
],
)

View File

@@ -0,0 +1,66 @@
# Set covering
An instance of set covering is composed of two entities: elements and sets; sets
cover a series of elements. The problem of set covering is about finding the
cheapest combination of sets that cover all the elements.
[More information on Wikipedia](https://en.wikipedia.org/wiki/Set_cover_problem).
* Solver: [`set_cover_heuristics.h`](set_cover_heuristics.h).
* Instance representation: [`set_cover_model.h`](set_cover_model.h).
* Instance parser: [`set_cover_reader.h`](set_cover_reader.h).
Create an instance:
```cpp
// If the elements are integers, a subset can be a std::vector<int> (in a pair
// along its cost).
std::vector<std::pair<std::vector<int>, int>> subsets = ...;
SetCoverModel model;
for (const auto [subset, subset_cost] : subsets) {
model.AddEmptySubset(subset_cost)
for (const int element : subset) {
model.AddElementToLastSubset(element);
}
}
SetCoverLedger ledger(&model);
```
Solve it using a MIP solver (guarantees to yield the optimum solution if it has
enough time to run):
```cpp
SetCoverMip mip(&ledger);
mip.SetTimeLimitInSeconds(10);
mip.NextSolution();
SubsetBoolVector best_choices = ledger.GetSolution();
LOG(INFO) << "Cost: " << ledger.cost();
```
A custom combination of heuristics (10,000 iterations of: clearing 10% of the
variables, running a Chvatal greedy descent, using steepest local search):
```cpp
Cost best_cost = std::numeric_limits<Cost>::max();
SubsetBoolVector best_choices = ledger.GetSolution();
for (int i = 0; i < 10000; ++i) {
ledger.LoadSolution(best_choices);
ClearRandomSubsets(0.1 * model.num_subsets().value(), &ledger);
GreedySolutionGenerator greedy(&ledger);
CHECK(greedy.NextSolution());
SteepestSearch steepest(&ledger);
CHECK(steepest.NextSolution(10000));
EXPECT_TRUE(ledger.CheckSolution());
if (ledger.cost() < best_cost) {
best_cost = ledger.cost();
best_choices = ledger.GetSolution();
LOG(INFO) << "Better cost: " << best_cost << " at iteration = " << i;
}
}
ledger.LoadSolution(best_choices);
LOG(INFO) << "Best cost: " << ledger.cost();
```

View File

@@ -14,9 +14,14 @@
#ifndef OR_TOOLS_SET_COVER_BASE_TYPES_H_
#define OR_TOOLS_SET_COVER_BASE_TYPES_H_
#include <cstddef>
#include <cstdint>
#include <type_traits>
#include <vector>
#include "absl/log/check.h"
#include "absl/time/time.h"
#include "absl/types/span.h"
#include "ortools/base/strong_int.h"
#include "ortools/base/strong_vector.h"
#include "ortools/base/timer.h"
@@ -66,6 +71,312 @@ using SparseRowView = util_intops::StrongVector<ElementIndex, SparseRow>;
using SubsetBoolVector = util_intops::StrongVector<SubsetIndex, bool>;
using ElementBoolVector = util_intops::StrongVector<ElementIndex, bool>;
// Maps from element to subset. Useful to compress the sparse row view.
using ElementToSubsetVector =
util_intops::StrongVector<ElementIndex, SubsetIndex>;
template <typename EntryIndex, typename Index>
class CompressedStrongVectorIterator;
// A compressed vector of strong indices (e.g. SubsetIndex, ElementIndex), with
// EntryIndex indicating the position in the vector (e.g. RowEntryIndex or
// ColumnEntryIndex).
// The vector is compressed in a variable-length encoding, where each element
// is encoded as a delta from the previous element.
// The vector is stored in a byte stream, which is addressed byte-by-byte, which
// is necessary to store variable-length integers.
// The variable-length encoding is little-endian base128 (LEB128) as used by
// protobufs. Since we only use non-negative integers, there is no need to
// encode the sign bit (so non-zigzag encoding or similar techniques).
//
// TODO(user):
// There is a lot of room for optimization in this class.
// - Use a bit-twiddling approach to encode and decode integers.
// - Use another simpler variable encoding (base on vu128) using a single 64-bit
// word and limited to storing 8 bytes with 7 bits per byte. A 64-bit word would
// contain 8*7 = 56 bits. This is enough for an index in an array with
// 2^56 = 7.2e16 elements, and more that the address space of current (2025)
// machines.
// - Access memory by 64-bit words instead of bytes.
// - Make Codecs a template parameter of CompressedStrongVector (?)
template <typename EntryIndex, typename Index>
class CompressedStrongVector {
public:
using iterator = CompressedStrongVectorIterator<EntryIndex, Index>;
using const_iterator = CompressedStrongVectorIterator<EntryIndex, Index>;
using value_type = Index;
explicit CompressedStrongVector() : memorized_value_(0), byte_stream_() {}
// Initializes the compressed vector from a strong vector.
// TODO(user): try to guess the size of the compressed vector and pre-allocate
// it. Experience shows it consumes between 1 and 2 bytes per Index on
// average.
explicit CompressedStrongVector(
const util_intops::StrongVector<EntryIndex, Index>& strong_vector)
: memorized_value_(0), byte_stream_() {
for (const Index x : strong_vector) {
push_back(x);
}
}
explicit CompressedStrongVector(absl::Span<const Index> vec)
: memorized_value_(0), byte_stream_() {
for (const Index x : vec) {
push_back(x);
}
}
// Appends an x to the vector in a compressed form.
void push_back(Index x) { EncodeInteger(x.value()); }
// Returns a reference to the underlying byte stream.
const std::vector<uint8_t>& byte_stream() const { return byte_stream_; }
// Returns the number of bytes needed to store the vector.
size_t size_in_bytes() const { return byte_stream_.size(); }
bool empty() const { return byte_stream_.empty(); }
// Reserves space for n bytes.
void Reserve(size_t n) { byte_stream_.reserve(n); }
// For range-for loops.
iterator begin() { return iterator(*this); }
iterator end() { return iterator(*this, true); }
// For const range-for loops.
const_iterator begin() const { return const_iterator(*this); }
const_iterator end() const { return const_iterator(*this, true); }
private:
void EncodeInteger(BaseInt x) {
BaseInt delta = x - memorized_value_; // Delta from previous value.
DCHECK_GE(delta, 0);
// Push the delta as a variable-length integer.
while (delta >= 128) {
// Keep the lower 7 bits of the delta and set the higher bit to 1 to mark
// the continuation of the variable-length encoding.
byte_stream_.push_back(static_cast<uint8_t>(delta | 0x80));
// Shift the delta by 7 bits to get the next value.
delta >>= 7;
}
// The final byte is less than 128, so its higher bit is zero, which marks
// the end of the variable-length encoding.
byte_stream_.push_back(static_cast<uint8_t>(delta));
// Do not forget to remember the last value.
memorized_value_ = x + kMinDelta;
}
// The last value memorized in the vector. It is defined as the last value
// appended to the vector plus a kMinDelta. It starts at zero.
BaseInt memorized_value_;
// The minimum difference between two consecutive elements in the vector.
static constexpr int64_t kMinDelta = 1;
// The byte stream.
std::vector<uint8_t> byte_stream_;
};
// Iterator for a compressed strong vector. There is no tool for decompressing
// a compressed strong vector, so this iterator is the only way to access the
// elements, always in order.
// The iterator is not thread-safe.
template <typename EntryIndex, typename Index>
class CompressedStrongVectorIterator {
public:
explicit CompressedStrongVectorIterator(
const CompressedStrongVector<EntryIndex, Index>& compressed_vector)
: compressed_vector_(compressed_vector),
memorized_value_(0),
pos_(0),
last_pos_(0) {}
explicit CompressedStrongVectorIterator(
const CompressedStrongVector<EntryIndex, Index>& compressed_vector,
bool at_end)
: compressed_vector_(compressed_vector),
memorized_value_(0),
pos_(0),
last_pos_(at_end ? compressed_vector.size_in_bytes() : 0) {}
bool at_end() const {
DCHECK_LE(last_pos_, compressed_vector_.size_in_bytes());
return last_pos_ == compressed_vector_.size_in_bytes();
}
Index operator*() { return DecodeInteger(); }
bool operator==(const CompressedStrongVectorIterator& other) const {
DCHECK_EQ(&compressed_vector_, &(other.compressed_vector_));
return last_pos_ == other.last_pos_;
}
bool operator!=(const CompressedStrongVectorIterator& other) const {
return !(*this == other);
}
CompressedStrongVectorIterator& operator++() {
last_pos_ = pos_;
return *this;
}
private:
Index DecodeInteger() {
// This can be made much faster by using a bit-twiddling approach.
const std::vector<uint8_t>& encoded = compressed_vector_.byte_stream();
BaseInt delta = 0;
int shift = 0;
for (pos_ = last_pos_, shift = 0;
encoded[pos_] >= 128 && pos_ < compressed_vector_.size_in_bytes();
shift += 7, ++pos_) {
delta |= static_cast<BaseInt>(encoded[pos_] & 0x7F) << shift;
}
delta |= static_cast<BaseInt>(encoded[pos_]) << shift;
++pos_;
memorized_value_ += Index(delta + kMinDelta);
return memorized_value_ - Index(kMinDelta);
}
// The compressed vector.
const CompressedStrongVector<EntryIndex, Index>& compressed_vector_;
// The last value memorized in the decoder. It is defined as the last value
// appended to the vector plus a kMinDelta. It starts at zero.
Index memorized_value_;
// The current position in the byte stream.
int64_t pos_;
// The last position in the byte stream.
int64_t last_pos_;
static constexpr int64_t kMinDelta = 1;
};
using CompressedColumn = CompressedStrongVector<ColumnEntryIndex, ElementIndex>;
using CompressedRow = CompressedStrongVector<RowEntryIndex, SubsetIndex>;
using CompressedColumnView =
util_intops::StrongVector<SubsetIndex, CompressedColumn>;
using CompressedRowView =
util_intops::StrongVector<ElementIndex, CompressedRow>;
using CompressedColumnIterator =
CompressedStrongVectorIterator<ColumnEntryIndex, ElementIndex>;
using CompressedRowIterator =
CompressedStrongVectorIterator<RowEntryIndex, SubsetIndex>;
template <typename Index>
class IndexRangeIterator;
// A range of indices, that can be iterated over. Useful if used in a range-for
// loop or as an IterableContainer.
template <typename Index>
class IndexRange {
public:
using iterator = IndexRangeIterator<Index>;
using const_iterator = IndexRangeIterator<Index>;
using value_type = Index;
IndexRange(Index start, Index end) : start_(start), end_(end) {}
explicit IndexRange(Index end) : start_(Index(0)), end_(end) {}
Index get_start() const { return start_; }
Index get_end() const { return end_; }
// For range-for loops.
iterator begin() { return iterator(*this); }
iterator end() { return iterator(*this, true); }
// For const range-for loops.
const_iterator begin() const { return const_iterator(*this); }
const_iterator end() const { return const_iterator(*this, true); }
private:
Index start_;
Index end_;
};
// Additional deduction guide
template <class Index>
IndexRange(Index a, Index b) -> IndexRange<Index>;
// The iterator for an IndexRange.
template <typename Index>
class IndexRangeIterator {
public:
explicit IndexRangeIterator(const IndexRange<Index>& range)
: range_(range), current_(range.get_start()) {}
IndexRangeIterator(const IndexRange<Index>& range, bool at_end)
: range_(range), current_(at_end ? range.get_end() : range.get_start()) {}
bool at_end() const { return current_ == range_.get_end(); }
Index operator*() { return current_; }
bool operator==(const IndexRangeIterator& other) const {
DCHECK_EQ(range_.get_start(), other.range_.get_start());
DCHECK_EQ(range_.get_end(), other.range_.get_end());
return current_ == other.current_;
}
bool operator!=(const IndexRangeIterator& other) const {
return !(*this == other);
}
IndexRangeIterator& operator++() {
++current_;
return *this;
}
private:
IndexRange<Index> range_;
Index current_;
};
// A container that can be iterated over, but does not own the data.
// The container can be either const or non-const.
// The container can be either a std::vector, a absl::Span, an IndexRange, a
// StrongVector or a CompressedStrongVector.
// We use the Curiously-Recurring Template Pattern (CRTP) to avoid having to
// specify the type of the container in the derived class, and to not lose
// runtime performance because of virtual functions.
template <typename T, typename Derived>
class IterableContainerBase {
public:
using value_type = typename std::decay_t<T>::value_type;
using iterator_type = decltype(std::begin(std::declval<T>()));
auto begin() { return static_cast<Derived*>(this)->begin_impl(); }
auto end() { return static_cast<Derived*>(this)->end_impl(); }
auto begin() const { return static_cast<const Derived*>(this)->begin_impl(); }
auto end() const { return static_cast<const Derived*>(this)->end_impl(); }
};
template <typename T>
class IterableContainer
: public IterableContainerBase<T, IterableContainer<T>> {
public:
explicit IterableContainer(const T& data_source) : data_(data_source) {}
auto begin_impl() { return data_.begin(); }
auto end_impl() { return data_.end(); }
auto begin_impl() const { return data_.begin(); }
auto end_impl() const { return data_.end(); }
private:
T data_;
};
// Additional deduction guide.
template <typename T>
IterableContainer(const T& data_source) -> IterableContainer<T>;
// Simple stopwatch class that enables recording the time spent in various
// functions in the code.
// It uses RAII to automatically record the time spent in the constructor and

View File

@@ -0,0 +1,200 @@
// Copyright 2010-2025 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/set_cover/base_types.h"
#include <cstddef>
#include <cstdint>
#include <random>
#include <vector>
#include "absl/types/span.h"
#include "benchmark/benchmark.h"
#include "gtest/gtest.h"
namespace operations_research {
namespace {
template <typename IndexType>
std::vector<IndexType> GenerateIncreasingVector(int start_value, int size) {
std::vector<IndexType> data(size);
std::iota(data.begin(), data.end(), IndexType(start_value));
return data;
}
template <typename ContainerType, typename IndexType>
class IterableContainerTest : public ::testing::Test {
public:
void SetUp() override {}
void TearDown() override {}
void CompareVectors(const ContainerType& data_source,
const std::vector<IndexType>& expected_data) {
IterableContainer<ContainerType> container(data_source);
std::vector<IndexType> actual_data;
for (const IndexType x : container) {
actual_data.push_back(x);
}
ASSERT_EQ(actual_data.size(), expected_data.size());
for (size_t i = 0; i < actual_data.size(); ++i) {
ASSERT_EQ(actual_data[i].value(), expected_data[i].value());
}
}
};
using VectorSubsetIndexTest =
IterableContainerTest<std::vector<SubsetIndex>, SubsetIndex>;
using VectorElementIndexTest =
IterableContainerTest<std::vector<ElementIndex>, ElementIndex>;
using SpanSubsetIndexTest =
IterableContainerTest<absl::Span<SubsetIndex>, SubsetIndex>;
using SpanElementIndexTest =
IterableContainerTest<absl::Span<ElementIndex>, ElementIndex>;
using IndexRangeSubsetIndexTest =
IterableContainerTest<IndexRange<SubsetIndex>, SubsetIndex>;
using IndexRangeElementIndexTest =
IterableContainerTest<IndexRange<ElementIndex>, ElementIndex>;
using SparseRowTest = IterableContainerTest<SparseRow, SubsetIndex>;
using SparseColumnTest = IterableContainerTest<SparseColumn, ElementIndex>;
using CompressedRowTest = IterableContainerTest<CompressedRow, SubsetIndex>;
using CompressedColumnTest =
IterableContainerTest<CompressedColumn, ElementIndex>;
TEST_F(VectorSubsetIndexTest, StdVectorSubsetIndex) {
std::vector<SubsetIndex> vec = {SubsetIndex(1), SubsetIndex(2),
SubsetIndex(3)};
std::vector<SubsetIndex> expected_data = vec;
CompareVectors(vec, expected_data);
}
TEST_F(VectorElementIndexTest, StdVectorElementIndex) {
std::vector<ElementIndex> vec = {ElementIndex(10), ElementIndex(20),
ElementIndex(30)};
std::vector<ElementIndex> expected_data = vec;
CompareVectors(vec, expected_data);
}
TEST_F(SpanSubsetIndexTest, AbslSpanSubsetIndex) {
std::vector<SubsetIndex> vec = {SubsetIndex(5), SubsetIndex(6),
SubsetIndex(7), SubsetIndex(8)};
absl::Span<SubsetIndex> span(vec);
CompareVectors(span, vec);
}
TEST_F(SpanElementIndexTest, AbslSpanElementIndex) {
std::vector<ElementIndex> vec = {ElementIndex(50), ElementIndex(60),
ElementIndex(70)};
absl::Span<ElementIndex> span(vec);
CompareVectors(span, vec);
}
TEST_F(IndexRangeSubsetIndexTest, IndexRangeSubsetIndex) {
IndexRange range(SubsetIndex(10), SubsetIndex(20));
std::vector<SubsetIndex> expected_data =
GenerateIncreasingVector<SubsetIndex>(10, 10);
CompareVectors(range, expected_data);
}
TEST_F(IndexRangeElementIndexTest, IndexRangeElementIndex) {
IndexRange range(ElementIndex(100), ElementIndex(110));
std::vector<ElementIndex> expected_data =
GenerateIncreasingVector<ElementIndex>(100, 10);
CompareVectors(range, expected_data);
}
TEST_F(SparseRowTest, StrongVectorSubsetIndex) {
SparseRow strong_vector{SubsetIndex(1), SubsetIndex(5), SubsetIndex(7)};
std::vector<SubsetIndex> expected_data = {SubsetIndex(1), SubsetIndex(5),
SubsetIndex(7)};
CompareVectors(strong_vector, expected_data);
}
TEST_F(SparseColumnTest, StrongVectorElementIndex) {
SparseColumn strong_vector{ElementIndex(10), ElementIndex(20),
ElementIndex(30)};
std::vector<ElementIndex> expected_data = {ElementIndex(10), ElementIndex(20),
ElementIndex(30)};
CompareVectors(strong_vector, expected_data);
}
TEST_F(CompressedRowTest, CompressedStrongVectorSubsetIndex) {
SparseRow original_vector{SubsetIndex(2), SubsetIndex(5), SubsetIndex(8),
SubsetIndex(12)};
CompressedRow compressed_vector(original_vector);
std::vector<SubsetIndex> expected_data = {SubsetIndex(2), SubsetIndex(5),
SubsetIndex(8), SubsetIndex(12)};
CompareVectors(compressed_vector, expected_data);
}
TEST_F(CompressedColumnTest, CompressedStrongVectorElementIndex) {
SparseColumn original_vector{ElementIndex(100), ElementIndex(105),
ElementIndex(111)};
CompressedStrongVector<ColumnEntryIndex, ElementIndex> compressed_vector(
original_vector);
std::vector<ElementIndex> expected_data = {
ElementIndex(100), ElementIndex(105), ElementIndex(111)};
CompareVectors(compressed_vector, expected_data);
}
SparseRow GenerateRandomSparseRow(size_t size, int64_t max_value) {
SparseRow sparse_row;
sparse_row.reserve(size);
std::mt19937_64 gen(std::random_device{}()); // Seed the generator
std::uniform_int_distribution<int64_t> dist(0, max_value);
SubsetIndex current_value(0);
for (size_t i = 0; i < size; ++i) {
current_value += SubsetIndex(dist(gen));
sparse_row.push_back(current_value);
}
return sparse_row;
}
static void BM_StrongVectorIteration(benchmark::State& state) {
const size_t size = state.range(0);
const int64_t delta_range = state.range(1);
SparseRow strong_vector = GenerateRandomSparseRow(size, delta_range);
for (auto _ : state) {
int64_t sum = 0;
for (const auto& x : strong_vector) {
sum += x.value();
}
benchmark::DoNotOptimize(sum); // Prevent optimization
}
}
BENCHMARK(BM_StrongVectorIteration)
->ArgsProduct({{100'000, 100'000'000}, {1 << 8, 1 << 16}});
static void BM_CompressedStrongVectorIteration(benchmark::State& state) {
const size_t size = state.range(0);
const int64_t delta_range = state.range(1);
CompressedRow compressed_vector(GenerateRandomSparseRow(size, delta_range));
for (auto _ : state) {
int64_t sum = 0;
for (const auto& x : compressed_vector) {
sum += x.value();
}
benchmark::DoNotOptimize(sum); // Prevent optimization
}
}
BENCHMARK(BM_CompressedStrongVectorIteration)
->ArgsProduct({
{100'000, 100'000'000},
{1 << 8, 1 << 16},
});
} // namespace
} // namespace operations_research

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,407 @@
// Copyright 2025 Francesco Cavaliere
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H
#define OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H
#include <absl/algorithm/container.h>
#include <absl/base/internal/pretty_function.h>
#include <absl/status/status.h>
#include <limits>
#include "ortools/set_cover/base_types.h"
#include "ortools/set_cover/set_cover_submodel.h"
#include "ortools/set_cover/set_cover_views.h"
namespace operations_research::scp {
// Implementation of:
// Caprara, Alberto, Matteo Fischetti, and Paolo Toth. 1999. “A Heuristic
// Method for the Set Covering Problem.” Operations Research 47 (5): 73043.
// https://www.jstor.org/stable/223097
//
// Hereafter referred to as CFT.
//
// SUMMARY
// The CFT algorithm is a heuristic approach to the set-covering problem. At its
// core, it combines a primal greedy heuristic with dual information obtained
// from the optimization of the Lagrangian relaxation of the problem.
// (Note: columns/subsets and rows/elements will be used interchangeably.)
//
// STRUCTURE
// The core of the algorithm is the 3Phase:
// 1. Subgradient optimization of the Lagrangian relaxation.
// 2. A primal greedy heuristic guided by the dual information.
// 3. Fixing some of the "best" columns (in terms of reduced costs) into the
// solution (diving).
// + Repeat until an exit criterion is met.
//
// The paper also considers an optional outer loop, which invokes the 3Phase
// process and then fixes some columns from the current best solution. This
// introduces two levels of diving: the outer loop fixes "primal-good" columns
// (based on the best solution), while the inner loop fixes "dual-good" columns
// (based on reduced costs).
// NOTE: The outer loop is not implemented in this version (yet - April 2025).
//
// Key characteristics of the algorithm:
//
// - The CFT algorithm is tailored for instances where the number of columns is
// significantly larger than the number of rows.
//
// - To improve efficiency, a core model approach is used. This involves
// selecting a small subset of columns based on their reduced costs, thereby
// substantially reducing the problem size handled in the internal steps.
//
// - Due to the use of the core model and column fixing, the algorithm rarely
// considers the entire problem. Instead, it operates on a small "window" of
// the problem. Efficiently managing this small window is a central aspect of
// any CFT implementation.
//
// - The core model scheme also enables an alternative implementation where the
// algorithm starts with a small model and progressively adds columns through
// a column-generation procedure. While this column generation procedure is
// problem-dependent and cannot be implemented here, the architecture of this
// implementation is designed to be extensible, allowing for such a procedure
// to be added in the future.
//
////////////////////////////////////////////////////////////////////////
////////////////////////// COMMON DEFINITIONS //////////////////////////
////////////////////////////////////////////////////////////////////////
// Small class to store the solution of a sub-model. It contains the cost and
// the subset list.
class Solution {
public:
Solution() = default;
Solution(const SubModel& model, const std::vector<SubsetIndex>& core_subsets);
double cost() const { return cost_; }
const std::vector<FullSubsetIndex>& subsets() const { return subsets_; }
void AddSubset(FullSubsetIndex subset, Cost cost) {
subsets_.push_back(subset);
cost_ += cost;
}
bool Empty() const { return subsets_.empty(); }
void Clear() {
cost_ = 0.0;
subsets_.clear();
}
private:
Cost cost_;
std::vector<FullSubsetIndex> subsets_;
};
// In the narrow scope of the CFT subgradient, there are often divisions
// between non-negative quantities (e.g., to compute a relative gap). In these
// specific cases, the denominator should always be greater than the
// numerator. This function checks that.
inline Cost DivideIfGE0(Cost numerator, Cost denominator) {
DCHECK_GE(numerator, -1e-6);
if (numerator < 1e-6) {
return 0.0;
}
return numerator / denominator;
}
// Dual information related to a SubModel.
// Stores multipliers, reduced costs, and the lower bound, and provides an
// interface that keeps them alligned.
class DualState {
public:
DualState() = default;
template <typename SubModelT>
DualState(const SubModelT& model)
: lower_bound_(.0),
multipliers_(model.num_elements(), .0),
reduced_costs_(model.subset_costs().begin(),
model.subset_costs().end()) {}
Cost lower_bound() const { return lower_bound_; }
const ElementCostVector& multipliers() const { return multipliers_; }
const SubsetCostVector& reduced_costs() const { return reduced_costs_; }
// NOTE: This function contains one of the two O(nnz) subgradient steps
template <typename SubModelT, typename Op>
void DualUpdate(const SubModelT& model, Op multiplier_operator) {
multipliers_.resize(model.num_elements());
reduced_costs_.resize(model.num_subsets());
lower_bound_ = .0;
// Update multipliers
for (ElementIndex i : model.ElementRange()) {
multiplier_operator(i, multipliers_[i]);
lower_bound_ += multipliers_[i];
DCHECK_GE(multipliers_[i], .0);
}
lower_bound_ += ComputeReducedCosts(model, multipliers_, reduced_costs_);
}
private:
// Single hot point to optimize once for the different use cases.
template <typename ModelT>
static Cost ComputeReducedCosts(const ModelT& model,
const ElementCostVector& multipliers,
SubsetCostVector& reduced_costs) {
// Compute new reduced costs (O(nnz))
Cost negative_sum = .0;
for (SubsetIndex j : model.SubsetRange()) {
reduced_costs[j] = model.subset_costs()[j];
for (ElementIndex i : model.columns()[j]) {
reduced_costs[j] -= multipliers[i];
}
if (reduced_costs[j] < .0) {
negative_sum += reduced_costs[j];
}
}
return negative_sum;
}
Cost lower_bound_;
ElementCostVector multipliers_;
SubsetCostVector reduced_costs_;
};
// Utility aggregate to store and pass around both primal and dual states.
struct PrimalDualState {
Solution solution;
DualState dual_state;
};
///////////////////////////////////////////////////////////////////////
///////////////////////////// SUBGRADIENT /////////////////////////////
///////////////////////////////////////////////////////////////////////
// Utilitiy aggregate used by the SubgradientOptimization procedure to
// communicate pass the needed information to the SubgradientCBs interface.
struct SubgradientContext {
const SubModel& model;
const DualState& current_dual_state;
const DualState& best_dual_state;
const Solution& best_solution;
const ElementCostVector& subgradient;
};
// Generic set of callbacks hooks used to specialized the behavior of the
// subgradient optimization
class SubgradientCBs {
public:
virtual bool ExitCondition(const SubgradientContext&) = 0;
virtual void RunHeuristic(const SubgradientContext&, Solution&) = 0;
virtual void ComputeMultipliersDelta(const SubgradientContext&,
ElementCostVector& delta_mults) = 0;
virtual bool UpdateCoreModel(SubModel&, PrimalDualState&) = 0;
virtual ~SubgradientCBs() = default;
};
// Subgradient optimization procedure. Optimizes the Lagrangian relaxation of
// the Set-Covering problem until a termination criterion si met.
void SubgradientOptimization(SubModel& core_model, SubgradientCBs& cbs,
PrimalDualState& best_state);
// Subgradient callbacks implementation focused on improving the current best
// dual bound.
class BoundCBs : public SubgradientCBs {
public:
static constexpr Cost kTol = 1e-6;
BoundCBs(const SubModel& model);
Cost step_size() const { return step_size_; }
bool ExitCondition(const SubgradientContext& context) override;
void ComputeMultipliersDelta(const SubgradientContext& context,
ElementCostVector& delta_mults) override;
void RunHeuristic(const SubgradientContext& context,
Solution& solution) override {}
bool UpdateCoreModel(SubModel& core_model,
PrimalDualState& best_state) override;
private:
void MakeMinimalCoverageSubgradient(const SubgradientContext& context,
ElementCostVector& subgradient);
private:
Cost squared_norm_;
// This addition implements a simplified stabilization technique inspired by
// works in the literature, such as:
// Frangioni, A., Gendron, B., Gorgone, E. (2015):
// "On the computational efficiency of subgradient methods: A case study in
// combinatorial optimization."
//
// The approach aims to reduce oscillations (zig-zagging) in the subgradient
// by using a "moving average" of the current and previous subgradients. The
// current subgradient is weighted by a factor of alpha, while the previous
// subgradients contribution is weighted by (1 - alpha). The parameter alpha
// is set to 0.5 by default but can be adjusted for tuning. The resulting
// stabilized subgradient is referred to as "direction" to distinguish it from
// the original subgradient.
Cost stabilization_coeff = 0.5; // Arbitrary from c4v4
ElementCostVector direction_;
ElementCostVector prev_direction_;
std::vector<SubsetIndex> lagrangian_solution_;
// Stopping condition
Cost prev_best_lb_;
BaseInt max_iter_countdown_;
BaseInt exit_test_countdown_;
BaseInt exit_test_period_;
BaseInt last_core_update_countdown_;
// Step size
void UpdateStepSize(Cost lower_bound);
Cost step_size_;
Cost last_min_lb_seen_;
Cost last_max_lb_seen_;
BaseInt step_size_update_countdown_;
BaseInt step_size_update_period_;
};
////////////////////////////////////////////////////////////////////////
/////////////////////// MULTIPLIERS BASED GREEDY ///////////////////////
////////////////////////////////////////////////////////////////////////
// Higher level function that return a function obtained by the dual multiplier
// based greedy.
Solution RunMultiplierBasedGreedy(
const SubModel& model, const DualState& dual_state,
Cost cost_cutoff = std::numeric_limits<BaseInt>::max());
// Lower level greedy function which creates or completes a "solution" (seen as
// set of subsets of the current SubModel) until a cutoff size or cost is
// reached.
// Note: since the cutoff might be reached, the returned solution might not be
// feasible.
Cost CoverGreedly(const SubModel& model, const DualState& dual_state,
Cost cost_cutoff, BaseInt size_cutoff,
std::vector<SubsetIndex>& sol_subsets);
///////////////////////////////////////////////////////////////////////
//////////////////////// THREE PHASE ALGORITHM ////////////////////////
///////////////////////////////////////////////////////////////////////
// Subgradient callbacks implementation focused on wandering near the optimal
// multipliers and invoke the multipliers based greedy heuristic at each
// iteration.
class HeuristicCBs : public SubgradientCBs {
public:
HeuristicCBs() : step_size_(0.1), countdown_(250) {};
void set_step_size(Cost step_size) { step_size_ = step_size; }
bool ExitCondition(const SubgradientContext& context) override {
Cost upper_bound =
context.best_solution.cost() - context.model.fixed_cost();
Cost lower_bound = context.best_dual_state.lower_bound();
return upper_bound - .999 < lower_bound || --countdown_ <= 0;
}
void RunHeuristic(const SubgradientContext& context,
Solution& solution) override;
void ComputeMultipliersDelta(const SubgradientContext& context,
ElementCostVector& delta_mults) override;
bool UpdateCoreModel(SubModel& model, PrimalDualState& state) override {
return false;
}
private:
Cost step_size_;
BaseInt countdown_;
};
PrimalDualState RunThreePhase(SubModel& model,
const Solution& init_solution = {});
///////////////////////////////////////////////////////////////////////
//////////////////////// FULL TO CORE PRICING /////////////////////////
///////////////////////////////////////////////////////////////////////
// CoreModel extractor. Stores a pointer to the full model and specilized
// UpdateCore in such a way to updated the SubModel (stored as base class) and
// focus the search on a small windows of the full model.
class FullToCoreModel : public SubModel {
using base = SubModel;
struct UpdateTrigger {
BaseInt countdown;
BaseInt period;
BaseInt max_period;
};
public:
FullToCoreModel(const Model* full_model);
Cost FixColumns(const std::vector<SubsetIndex>& columns_to_fix) override;
bool UpdateCore(PrimalDualState& core_state) override;
void ResetPricingPeriod();
const DualState& best_dual_state() const { return best_dual_state_; }
private:
decltype(auto) IsFocusCol(FullSubsetIndex j) {
return is_focus_col_[static_cast<SubsetIndex>(j)];
}
decltype(auto) IsFocusRow(FullElementIndex i) {
return is_focus_row_[static_cast<ElementIndex>(i)];
}
void UpdatePricingPeriod(const DualState& full_dual_state,
const PrimalDualState& core_state);
void UpdateDualState(const DualState& core_dual_state,
DualState& full_dual_state, DualState& best_dual_state);
// Views are not composable (for now), so we can either access the full_model
// with the strongly typed view or with the filtered view.
// Access the full model filtered by the current columns fixed.
FilterModelView FixingFullModelView() const {
return FilterModelView(full_model_, &is_focus_col_, &is_focus_row_,
num_subsets_, num_elements_);
}
// Access the full model with the strongly typed view.
StrongModelView StrongTypedFullModelView() const {
return StrongModelView(full_model_);
}
bool FullToSubModelInvariantCheck();
private:
const Model* full_model_;
// Note: The `is_focus_col_` vector duplicates information already present in
// `SubModelView::cols_sizes_`. However, it does not overlap with any data
// stored in `CoreModel`. Since `CoreModel` is expected to be the primary use
// case, this vector is explicitly maintained here to ensure compatibility.
SubsetBoolVector is_focus_col_;
// Note: The `is_focus_row_` vector is functionally redundant with either
// `CoreModel::full2core_row_map_` or `SubModelView::rows_sizes_`. These
// existing structures could be used to create the filtered view of the full
// model. However, doing so would require generalizing the current view
// system to work with generic functors instead of vectors of integral types.
// Since the number of elements is assumed to be not prohibitive, a simpler
// implementation that avoids this memory optimization was preferred.
ElementBoolVector is_focus_row_;
BaseInt num_subsets_;
BaseInt num_elements_;
DualState full_dual_state_;
DualState best_dual_state_;
BaseInt update_countdown_;
BaseInt update_period_;
BaseInt update_max_period_;
};
// Coverage counter to decide the number of columns to keep in the core model.
static constexpr BaseInt kMinCov = 5;
} // namespace operations_research::scp
#endif /* OR_TOOLS_ORTOOLS_SET_COVER_SET_COVER_CFT_H */

View File

@@ -42,6 +42,10 @@ static constexpr double kInfinity = std::numeric_limits<float>::infinity();
using CL = SetCoverInvariant::ConsistencyLevel;
bool SetCoverSolutionGenerator::CheckInvariantConsistency() const {
return inv_->CheckConsistency(consistency_level_);
}
// TrivialSolutionGenerator.
bool TrivialSolutionGenerator::NextSolution(
@@ -505,10 +509,8 @@ bool SteepestSearch::NextSolution(const SubsetBoolVector& in_focus) {
DCHECK(inv()->ComputeIsRedundant(best_subset));
DCHECK_GT(costs[best_subset], 0.0);
inv()->Deselect(best_subset, CL::kFreeAndUncovered);
for (IntersectingSubsetsIterator it(*model(), best_subset); !it.at_end();
++it) {
const SubsetIndex subset = *it;
for (const SubsetIndex subset :
IntersectingSubsetsRange(*model(), best_subset)) {
if (!inv()->ComputeIsRedundant(subset)) {
pq.Remove(subset.value());
}
@@ -544,7 +546,7 @@ bool LazySteepestSearch::NextSolution(const SubsetBoolVector& in_focus) {
}
std::vector<BaseInt> cost_permutation(selected_costs.size());
std::iota(cost_permutation.begin(), cost_permutation.end(), 0);
// TODO(user): use RadixSort and sort on Cost,we need doubles and payloads.
// TODO(user): use RadixSort with doubles and payloads.
std::sort(cost_permutation.begin(), cost_permutation.end(),
[&selected_costs](const BaseInt i, const BaseInt j) {
if (selected_costs[i] == selected_costs[j]) {
@@ -826,29 +828,35 @@ std::vector<SubsetIndex> ClearRandomSubsets(BaseInt num_subsets,
}
std::vector<SubsetIndex> ClearRandomSubsets(absl::Span<const SubsetIndex> focus,
BaseInt num_subsets,
BaseInt num_subsets_to_choose,
SetCoverInvariant* inv) {
num_subsets = std::min<BaseInt>(num_subsets, focus.size());
CHECK_GE(num_subsets, 0);
num_subsets_to_choose =
std::min<BaseInt>(num_subsets_to_choose, focus.size());
CHECK_GE(num_subsets_to_choose, 0);
std::vector<SubsetIndex> chosen_indices;
for (const SubsetIndex subset : focus) {
if (inv->is_selected()[subset]) {
chosen_indices.push_back(subset);
}
}
SampleSubsets(&chosen_indices, num_subsets);
SampleSubsets(&chosen_indices, num_subsets_to_choose);
BaseInt num_deselected = 0;
for (const SubsetIndex subset : chosen_indices) {
inv->Deselect(subset, CL::kCostAndCoverage);
++num_deselected;
for (IntersectingSubsetsIterator it(*inv->model(), subset); !it.at_end();
++it) {
if (!inv->is_selected()[subset]) continue;
if (inv->is_selected()[subset]) { // subset may have been deselected in a
// previous iteration.
inv->Deselect(subset, CL::kCostAndCoverage);
++num_deselected;
}
// Note that num_deselected may exceed num_subsets by more than 1.
if (num_deselected > num_subsets) break;
for (const SubsetIndex connected_subset :
IntersectingSubsetsRange(*inv->model(), subset)) {
// connected_subset may have been deselected in a previous iteration.
if (inv->is_selected()[connected_subset]) {
inv->Deselect(connected_subset, CL::kCostAndCoverage);
++num_deselected;
}
}
// Note that num_deselected may exceed num_subsets_to_choose by more than 1.
if (num_deselected > num_subsets_to_choose) break;
}
return chosen_indices;
}

View File

@@ -14,8 +14,6 @@
#ifndef OR_TOOLS_SET_COVER_SET_COVER_HEURISTICS_H_
#define OR_TOOLS_SET_COVER_SET_COVER_HEURISTICS_H_
#include <stdbool.h>
#include <cstdint>
#include <limits>
#include <string>
@@ -55,10 +53,12 @@ class SetCoverSolutionGenerator {
// By default, the maximum number of iterations is set to infinity, and the
// maximum time in seconds is set to infinity as well (and the time limit is
// not yet implemented).
SetCoverSolutionGenerator(SetCoverInvariant* inv,
absl::string_view class_name,
absl::string_view name)
SetCoverSolutionGenerator(
SetCoverInvariant* inv,
SetCoverInvariant::ConsistencyLevel consistency_level,
absl::string_view class_name, absl::string_view name)
: run_time_(absl::ZeroDuration()),
consistency_level_(consistency_level),
inv_(inv),
class_name_(class_name),
name_(name),
@@ -126,6 +126,8 @@ class SetCoverSolutionGenerator {
// Same as above, but with a vector of Booleans as focus.
virtual bool NextSolution(const SubsetBoolVector& in_focus) = 0;
bool CheckInvariantConsistency() const;
protected:
// Accessors.
SetCoverModel* model() const { return inv_->model(); }
@@ -137,6 +139,9 @@ class SetCoverSolutionGenerator {
// run_time_ is an abstract duration for the time spent in NextSolution().
absl::Duration run_time_;
// The consistency needed by the solution generator.
SetCoverInvariant::ConsistencyLevel consistency_level_;
private:
// The data structure that will maintain the invariant for the model.
SetCoverInvariant* inv_;
@@ -165,10 +170,11 @@ class SetCoverSolutionGenerator {
// subset indices if needed.
class SubsetListBasedSolutionGenerator : public SetCoverSolutionGenerator {
public:
explicit SubsetListBasedSolutionGenerator(SetCoverInvariant* inv,
absl::string_view class_name,
absl::string_view name)
: SetCoverSolutionGenerator(inv, class_name, name) {}
explicit SubsetListBasedSolutionGenerator(
SetCoverInvariant* inv,
SetCoverInvariant::ConsistencyLevel consistency_level,
absl::string_view class_name, absl::string_view name)
: SetCoverSolutionGenerator(inv, consistency_level, class_name, name) {}
bool NextSolution(absl::Span<const SubsetIndex> _) override { return false; }
@@ -201,10 +207,11 @@ class SubsetListBasedSolutionGenerator : public SetCoverSolutionGenerator {
// Booleans if needed.
class BoolVectorBasedSolutionGenerator : public SetCoverSolutionGenerator {
public:
explicit BoolVectorBasedSolutionGenerator(SetCoverInvariant* inv,
absl::string_view class_name,
absl::string_view name)
: SetCoverSolutionGenerator(inv, class_name, name) {}
explicit BoolVectorBasedSolutionGenerator(
SetCoverInvariant* inv,
SetCoverInvariant::ConsistencyLevel consistency_level,
absl::string_view class_name, absl::string_view name)
: SetCoverSolutionGenerator(inv, consistency_level, class_name, name) {}
bool NextSolution(const SubsetBoolVector& _) override { return false; }
@@ -241,7 +248,9 @@ class TrivialSolutionGenerator : public SubsetListBasedSolutionGenerator {
: TrivialSolutionGenerator(inv, "TrivialGenerator") {}
TrivialSolutionGenerator(SetCoverInvariant* inv, absl::string_view name)
: SubsetListBasedSolutionGenerator(inv, "TrivialGenerator", name) {}
: SubsetListBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kFreeAndUncovered,
"TrivialGenerator", name) {}
using SubsetListBasedSolutionGenerator::NextSolution;
bool NextSolution(absl::Span<const SubsetIndex> focus) final;
@@ -260,7 +269,9 @@ class RandomSolutionGenerator : public SubsetListBasedSolutionGenerator {
: RandomSolutionGenerator(inv, "RandomGenerator") {}
RandomSolutionGenerator(SetCoverInvariant* inv, absl::string_view name)
: SubsetListBasedSolutionGenerator(inv, "RandomGenerator", name) {}
: SubsetListBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kFreeAndUncovered,
"RandomGenerator", name) {}
using SubsetListBasedSolutionGenerator::NextSolution;
bool NextSolution(absl::Span<const SubsetIndex> focus) final;
@@ -293,7 +304,9 @@ class GreedySolutionGenerator : public SubsetListBasedSolutionGenerator {
: GreedySolutionGenerator(inv, "GreedyGenerator") {}
GreedySolutionGenerator(SetCoverInvariant* inv, absl::string_view name)
: SubsetListBasedSolutionGenerator(inv, "GreedyGenerator", name) {}
: SubsetListBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kFreeAndUncovered,
"GreedyGenerator", name) {}
using SubsetListBasedSolutionGenerator::NextSolution;
bool NextSolution(absl::Span<const SubsetIndex> focus) final;
@@ -314,7 +327,9 @@ class ElementDegreeSolutionGenerator : public BoolVectorBasedSolutionGenerator {
: ElementDegreeSolutionGenerator(inv, "ElementDegreeGenerator") {}
ElementDegreeSolutionGenerator(SetCoverInvariant* inv, absl::string_view name)
: BoolVectorBasedSolutionGenerator(inv, "ElementDegreeGenerator", name) {}
: BoolVectorBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kFreeAndUncovered,
"ElementDegreeGenerator", name) {}
using BoolVectorBasedSolutionGenerator::NextSolution;
bool NextSolution(const SubsetBoolVector& in_focus) final;
@@ -337,8 +352,9 @@ class LazyElementDegreeSolutionGenerator
LazyElementDegreeSolutionGenerator(SetCoverInvariant* inv,
absl::string_view name)
: BoolVectorBasedSolutionGenerator(inv, "LazyElementDegreeGenerator",
name) {}
: BoolVectorBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kCostAndCoverage,
"LazyElementDegreeGenerator", name) {}
using BoolVectorBasedSolutionGenerator::NextSolution;
bool NextSolution(const SubsetBoolVector& in_focus) final;
@@ -358,7 +374,9 @@ class SteepestSearch : public BoolVectorBasedSolutionGenerator {
: SteepestSearch(inv, "SteepestSearch") {}
SteepestSearch(SetCoverInvariant* inv, absl::string_view name)
: BoolVectorBasedSolutionGenerator(inv, "SteepestSearch", name) {}
: BoolVectorBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kFreeAndUncovered,
"SteepestSearch", name) {}
using BoolVectorBasedSolutionGenerator::NextSolution;
bool NextSolution(const SubsetBoolVector& in_focus) final;
@@ -375,7 +393,9 @@ class LazySteepestSearch : public BoolVectorBasedSolutionGenerator {
: LazySteepestSearch(inv, "LazySteepestSearch") {}
LazySteepestSearch(SetCoverInvariant* inv, absl::string_view name)
: BoolVectorBasedSolutionGenerator(inv, "LazySteepestSearch", name) {}
: BoolVectorBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kCostAndCoverage,
"LazySteepestSearch", name) {}
using BoolVectorBasedSolutionGenerator::NextSolution;
bool NextSolution(const SubsetBoolVector& in_focus) final;
@@ -459,7 +479,9 @@ class GuidedTabuSearch : public SubsetListBasedSolutionGenerator {
: GuidedTabuSearch(inv, "GuidedTabuSearch") {}
GuidedTabuSearch(SetCoverInvariant* inv, absl::string_view name)
: SubsetListBasedSolutionGenerator(inv, "GuidedTabuSearch", name),
: SubsetListBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kFreeAndUncovered,
"GuidedTabuSearch", name),
lagrangian_factor_(kDefaultLagrangianFactor),
penalty_factor_(kDefaultPenaltyFactor),
epsilon_(kDefaultEpsilon),
@@ -549,7 +571,9 @@ class GuidedLocalSearch : public SubsetListBasedSolutionGenerator {
: GuidedLocalSearch(inv, "GuidedLocalSearch") {}
GuidedLocalSearch(SetCoverInvariant* inv, absl::string_view name)
: SubsetListBasedSolutionGenerator(inv, "GuidedLocalSearch", name),
: SubsetListBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kRedundancy,
"GuidedLocalSearch", name),
epsilon_(kDefaultEpsilon),
alpha_(kDefaultAlpha) {
Initialize();

View File

@@ -341,15 +341,17 @@ BaseInt SetCoverInvariant::ComputeNumFreeElements(SubsetIndex subset) const {
return num_free_elements;
}
void SetCoverInvariant::Select(SubsetIndex subset,
bool SetCoverInvariant::Select(SubsetIndex subset,
ConsistencyLevel target_consistency) {
if (is_selected_[subset]) return false;
const bool update_redundancy_info = target_consistency >= CL::kRedundancy;
if (update_redundancy_info) {
ClearRemovabilityInformation();
}
consistency_level_ = std::min(consistency_level_, target_consistency);
DVLOG(1) << "Selecting subset " << subset;
DCHECK(!is_selected_[subset]);
DCHECK(CheckConsistency(target_consistency));
trace_.push_back(SetCoverDecision(subset, true));
is_selected_[subset] = true;
@@ -362,7 +364,7 @@ void SetCoverInvariant::Select(SubsetIndex subset,
for (const ElementIndex element : columns[subset]) {
++coverage_[element];
}
return;
return true;
}
for (const ElementIndex element : columns[subset]) {
if (coverage_[element] == 0) {
@@ -402,10 +404,14 @@ void SetCoverInvariant::Select(SubsetIndex subset,
}
DCHECK_EQ(num_free_elements_[subset], 0);
DCHECK(CheckConsistency(target_consistency));
return true;
}
void SetCoverInvariant::Deselect(SubsetIndex subset,
bool SetCoverInvariant::Deselect(SubsetIndex subset,
ConsistencyLevel target_consistency) {
// NOMUTANTS -- This is a short-circuit to avoid doing any work
// when the subset is already deselected.
if (!is_selected_[subset]) return false;
DCHECK(CheckConsistency(target_consistency));
const bool update_redundancy_info = target_consistency >= CL::kRedundancy;
if (update_redundancy_info) {
@@ -426,7 +432,7 @@ void SetCoverInvariant::Deselect(SubsetIndex subset,
for (const ElementIndex element : columns[subset]) {
--coverage_[element];
}
return;
return true;
}
// This is a dissymmetry with Select, and only maintained in
// consistency level kFreeAndUncovered and above.
@@ -461,6 +467,7 @@ void SetCoverInvariant::Deselect(SubsetIndex subset,
// nor meaning in adding it a list of removable or non-removable
// subsets. This is a dissymmetry with Select.
DCHECK(CheckConsistency(target_consistency));
return true;
}
SetCoverSolutionResponse SetCoverInvariant::ExportSolutionAsProto() const {

View File

@@ -227,11 +227,18 @@ class SetCoverInvariant {
// Includes subset in the solution by setting is_selected_[subset] to true
// and incrementally updating the invariant to the given consistency level.
void Select(SubsetIndex subset, ConsistencyLevel consistency);
// Returns false if the subset is already selected.
// This allows a programming style where Select is embedded in a DCHECK, or
// to write `if (Select(subset, consistency)) { ... }`, which is more readable
// than `if (!is_selected_[subset]) { Select(subset, consistency); ... }`
// The above is a good programming style for example to count how many
// elements were actually set.
bool Select(SubsetIndex subset, ConsistencyLevel consistency);
// Excludes subset from the solution by setting is_selected_[subset] to false
// and incrementally updating the invariant to the given consistency level.
void Deselect(SubsetIndex subset, ConsistencyLevel consistency);
// Returns false if the subset is already deselected.
bool Deselect(SubsetIndex subset, ConsistencyLevel consistency);
// Returns the current solution as a proto.
SetCoverSolutionResponse ExportSolutionAsProto() const;

View File

@@ -16,7 +16,6 @@
#include <memory>
#include <tuple>
#include <vector>
#include "absl/strings/string_view.h"
#include "absl/types/span.h"
@@ -24,7 +23,6 @@
#include "ortools/set_cover/base_types.h"
#include "ortools/set_cover/set_cover_heuristics.h"
#include "ortools/set_cover/set_cover_invariant.h"
#include "ortools/set_cover/set_cover_model.h"
namespace operations_research {
@@ -55,7 +53,9 @@ class SetCoverLagrangian : public SubsetListBasedSolutionGenerator {
: SetCoverLagrangian(inv, "Lagrangian") {}
SetCoverLagrangian(SetCoverInvariant* inv, const absl::string_view name)
: SubsetListBasedSolutionGenerator(inv, "Lagrangian", name),
: SubsetListBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kInconsistent,
"Lagrangian", name),
num_threads_(1),
thread_pool_(nullptr) {}

View File

@@ -135,7 +135,8 @@ bool SetCoverMip::NextSolution(absl::Span<const SubsetIndex> focus) {
if (use_integers_) {
using CL = SetCoverInvariant::ConsistencyLevel;
for (const SubsetIndex subset : focus) {
if (vars[subset]->solution_value() > 0.9) {
if (vars[subset]->solution_value() > 0.9 &&
!inv()->is_selected()[subset]) {
inv()->Select(subset, CL::kCostAndCoverage);
}
}

View File

@@ -37,7 +37,9 @@ class SetCoverMip : public SubsetListBasedSolutionGenerator {
: SetCoverMip(inv, "SetCoverMip") {}
SetCoverMip(SetCoverInvariant* inv, absl::string_view name)
: SubsetListBasedSolutionGenerator(inv, "Mip", name),
: SubsetListBasedSolutionGenerator(
inv, SetCoverInvariant::ConsistencyLevel::kCostAndCoverage, "Mip",
name),
mip_solver_(SetCoverMipSolver::SCIP),
use_integers_(true) {}
@@ -48,6 +50,9 @@ class SetCoverMip : public SubsetListBasedSolutionGenerator {
SetCoverMip& UseIntegers(bool use_integers) {
use_integers_ = use_integers;
consistency_level_ =
use_integers_ ? SetCoverInvariant::ConsistencyLevel::kCostAndCoverage
: SetCoverInvariant::ConsistencyLevel::kInconsistent;
return *this;
}

View File

@@ -20,6 +20,7 @@
#include "absl/log/check.h"
#include "absl/strings/string_view.h"
#include "absl/time/time.h"
#include "absl/types/span.h"
#include "ortools/base/strong_int.h"
#include "ortools/base/strong_vector.h"
@@ -395,7 +396,7 @@ class SetCoverModel {
// Vector of rows. Each row corresponds to an element and contains the
// subsets containing the element.
// The size is exactly the same as for columns_.
// The size is exactly the same as for columns_ (or there would be a bug.)
SparseRowView rows_;
// Vector of indices from 0 to columns.size() - 1. (Like std::iota, but built
@@ -407,52 +408,78 @@ class SetCoverModel {
std::vector<SubsetIndex> all_subsets_;
};
// The IntersectingSubsetsIterator is a forward iterator that returns the next
// IntersectingSubsetsIterator is a forward iterator that returns the next
// intersecting subset for a fixed seed_subset.
// The iterator is initialized with a model and a seed_subset and
// allows a speedup in getting the intersecting subsets
// by not storing them in memory.
// The iterator is at the end when the last intersecting subset has been
// returned.
// TODO(user): Add the possibility for range-for loops.
class IntersectingSubsetsIterator {
public:
IntersectingSubsetsIterator(const SetCoverModel& model,
SubsetIndex seed_subset)
: intersecting_subset_(-1),
element_entry_(0),
subset_entry_(0),
SubsetIndex seed_subset, bool at_end = false)
: model_(model),
seed_subset_(seed_subset),
model_(model),
subset_seen_(model_.columns().size(), false) {
seed_column_(model_.columns()[seed_subset]),
seed_column_size_(ColumnEntryIndex(seed_column_.size())),
intersecting_subset_(0), // meaningless, will be overwritten.
element_entry_(0),
rows_(model_.rows()),
subset_seen_() {
// For iterator to be as light as possible when created, we do not reserve
// space for the subset_seen_ vector, and we do not initialize it. This
// is done to avoid the overhead of creating the vector and initializing it
// when the iterator is created. The vector is created on the first call to
// the ++ operator.
DCHECK(model_.row_view_is_valid());
subset_seen_[seed_subset] = true; // Avoid iterating on `seed_subset`.
++(*this); // Move to the first intersecting subset.
if (at_end) {
element_entry_ = seed_column_size_;
return;
}
for (; element_entry_ < seed_column_size_; ++element_entry_) {
const ElementIndex current_element = seed_column_[element_entry_];
const SparseRow& current_row = rows_[current_element];
const RowEntryIndex current_row_size = RowEntryIndex(current_row.size());
for (; subset_entry_ < current_row_size; ++subset_entry_) {
if (intersecting_subset_ == seed_subset_) continue;
intersecting_subset_ = current_row[subset_entry_];
return;
}
subset_entry_ = RowEntryIndex(0); // 'carriage-return'
}
}
// Returns (true) whether the iterator is at the end.
bool at_end() const {
return element_entry_.value() == model_.columns()[seed_subset_].size();
}
bool at_end() const { return element_entry_ == seed_column_size_; }
// Returns the intersecting subset.
SubsetIndex operator*() const { return intersecting_subset_; }
// Move the iterator to the next intersecting subset.
// Disequality operator.
bool operator!=(const IntersectingSubsetsIterator& other) const {
return element_entry_ != other.element_entry_ ||
subset_entry_ != other.subset_entry_ ||
seed_subset_ != other.seed_subset_;
}
// Advances the iterator to the next intersecting subset.
IntersectingSubsetsIterator& operator++() {
DCHECK(model_.row_view_is_valid());
DCHECK(!at_end());
const SparseRowView& rows = model_.rows();
const SparseColumn& column = model_.columns()[seed_subset_];
const ColumnEntryIndex column_size = ColumnEntryIndex(column.size());
for (; element_entry_ < column_size; ++element_entry_) {
const ElementIndex current_element = column[element_entry_];
const SparseRow& current_row = rows[current_element];
DCHECK(!at_end()) << "element_entry_ = " << element_entry_
<< " subset_entry_ = " << subset_entry_
<< " seed_column_size_ = " << seed_column_size_;
if (subset_seen_.empty()) {
subset_seen_.resize(model_.num_subsets(), false);
subset_seen_[seed_subset_] = true;
}
subset_seen_[intersecting_subset_] = true;
for (; element_entry_ < seed_column_size_; ++element_entry_) {
const ElementIndex current_element = seed_column_[element_entry_];
const SparseRow& current_row = rows_[current_element];
const RowEntryIndex current_row_size = RowEntryIndex(current_row.size());
for (; subset_entry_ < current_row_size; ++subset_entry_) {
intersecting_subset_ = current_row[subset_entry_];
if (!subset_seen_[intersecting_subset_]) {
subset_seen_[intersecting_subset_] = true;
return *this;
}
}
@@ -462,6 +489,18 @@ class IntersectingSubsetsIterator {
}
private:
// The model to which the iterator is applying.
const SetCoverModel& model_;
// The seed subset.
SubsetIndex seed_subset_;
// A reference to the column of the seed subset, kept here for ease of access.
const SparseColumn& seed_column_;
// The size of the column of the seed subset.
ColumnEntryIndex seed_column_size_;
// The intersecting subset.
SubsetIndex intersecting_subset_;
@@ -471,16 +510,46 @@ class IntersectingSubsetsIterator {
// The position of the entry in the row corresponding to `element_entry`.
RowEntryIndex subset_entry_;
// The seed subset.
SubsetIndex seed_subset_;
// The model to which the iterator is applying.
const SetCoverModel& model_;
// A reference to the rows of the model, kept here for ease of access.
const SparseRowView& rows_;
// A vector of Booleans indicating whether the current subset has been
// already seen by the iterator.
SubsetBoolVector subset_seen_;
};
// IntersectingSubsetsRange is a range of intersecting subsets for a fixed seed
// subset. Can be used with range-based for-loops.
class IntersectingSubsetsRange {
public:
using iterator = IntersectingSubsetsIterator;
using const_iterator = IntersectingSubsetsIterator;
IntersectingSubsetsRange(const SetCoverModel& model, SubsetIndex seed_subset)
: model_(model), seed_subset_(seed_subset) {}
iterator begin() { return IntersectingSubsetsIterator(model_, seed_subset_); }
iterator end() {
return IntersectingSubsetsIterator(model_, seed_subset_, true);
}
const_iterator begin() const {
return IntersectingSubsetsIterator(model_, seed_subset_);
}
const_iterator end() const {
return IntersectingSubsetsIterator(model_, seed_subset_, true);
}
private:
// The model to which the range is applying.
const SetCoverModel& model_;
// The seed subset for which the intersecting subsets are computed.
SubsetIndex seed_subset_;
};
} // namespace operations_research
#endif // OR_TOOLS_SET_COVER_SET_COVER_MODEL_H_

View File

@@ -95,8 +95,9 @@ ABSL_FLAG(int, max_elements_for_classic, 5000,
ABSL_FLAG(bool, benchmarks, false, "Run benchmarks.");
ABSL_FLAG(std::string, benchmarks_dir, "", "Benchmarks directory.");
ABSL_FLAG(bool, collate_scp, false, "Collate the SCP benchmarks.");
// TODO(user): Add a flags to:
// TODO(user): Add flags to:
// - Choose problems by name or by size: filter_name, max_elements, max_subsets.
// - Exclude problems by name: exclude_name.
// - Choose which solution generators to run.
@@ -124,18 +125,14 @@ void LogStats(const SetCoverModel& model) {
LOG(INFO) << header << model.ToVerboseString(sep);
LOG(INFO) << header << "cost" << sep
<< model.ComputeCostStats().ToVerboseString(sep);
LOG(INFO) << header << "row stats" << sep
LOG(INFO) << header << "row size stats" << sep
<< model.ComputeRowStats().ToVerboseString(sep);
LOG(INFO) << header << "row size deciles" << sep
<< absl::StrJoin(model.ComputeRowDeciles(), sep);
LOG(INFO) << header << "row delta byte size stats" << sep
<< model.ComputeRowDeltaSizeStats().ToVerboseString(sep);
LOG(INFO) << header << "column stats" << sep
LOG(INFO) << header << "column size stats" << sep
<< model.ComputeColumnStats().ToVerboseString(sep);
LOG(INFO) << header << "column size deciles" << sep
<< absl::StrJoin(model.ComputeColumnDeciles(), sep);
LOG(INFO) << header << "column delta byte size stats" << sep
<< model.ComputeColumnDeltaSizeStats().ToVerboseString(sep);
LOG(INFO) << header << "num_singleton_rows" << sep
<< model.ComputeNumSingletonRows() << sep << "num_singleton_columns"
<< sep << model.ComputeNumSingletonColumns();
@@ -391,10 +388,11 @@ std::string FormatModelAndStats(const SetCoverModel& model) {
} else { // CSV
const std::string header =
absl::StrCat(Separator(), model.name(), Separator());
return absl::StrCat(
header, model.ToString(Separator()), Eol(), header, "column stats",
Separator(), FormatStats(model.ComputeColumnStats()), Eol(), header,
"row stats", Separator(), FormatStats(model.ComputeRowStats()));
return absl::StrCat(header, model.ToString(Separator()), Eol(), header,
"column size stats", Separator(),
FormatStats(model.ComputeColumnStats()), Eol(), header,
"row size stats", Separator(),
FormatStats(model.ComputeRowStats()), Eol());
}
}
@@ -502,17 +500,36 @@ std::vector<BenchmarksTableRow> BenchmarksTable() {
// directory specified by the --benchmarks_dir flag.
// Use a macro to be able to compute the size of the array at compile time.
#define BUILD_VECTOR(files) BuildVector(files, ABSL_ARRAYSIZE(files))
return std::vector<BenchmarksTableRow>{
{"scp-orlib", BUILD_VECTOR(kScp4To6Files), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpAToEFiles), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpNrFiles), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpClrFiles), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpCycFiles), FileFormat::ORLIB},
{"scp-rail", BUILD_VECTOR(kRailFiles), FileFormat::RAIL},
{"scp-wedelin", BUILD_VECTOR(kWedelinFiles), FileFormat::ORLIB},
{"scp-balas", BUILD_VECTOR(kBalasFiles), FileFormat::ORLIB},
{"scp-fimi", BUILD_VECTOR(kFimiFiles), FileFormat::FIMI},
};
std::vector<BenchmarksTableRow> result;
std::vector<std::string> all_scp_files;
if (absl::GetFlag(FLAGS_collate_scp)) {
all_scp_files = BUILD_VECTOR(kScp4To6Files);
all_scp_files.insert(all_scp_files.end(), kScpAToEFiles,
kScpAToEFiles + ABSL_ARRAYSIZE(kScpAToEFiles));
all_scp_files.insert(all_scp_files.end(), kScpNrFiles,
kScpNrFiles + ABSL_ARRAYSIZE(kScpNrFiles));
all_scp_files.insert(all_scp_files.end(), kScpCycFiles,
kScpCycFiles + ABSL_ARRAYSIZE(kScpCycFiles));
result = {{"scp-orlib", all_scp_files, FileFormat::ORLIB}};
} else {
result = {
{"scp-orlib", BUILD_VECTOR(kScp4To6Files), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpAToEFiles), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpNrFiles), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpClrFiles), FileFormat::ORLIB},
{"scp-orlib", BUILD_VECTOR(kScpCycFiles), FileFormat::ORLIB},
};
}
result.insert(
result.end(),
{
{"scp-rail", BUILD_VECTOR(kRailFiles), FileFormat::RAIL},
{"scp-wedelin", BUILD_VECTOR(kWedelinFiles), FileFormat::ORLIB},
{"scp-balas", BUILD_VECTOR(kBalasFiles), FileFormat::ORLIB},
{"scp-fimi", BUILD_VECTOR(kFimiFiles), FileFormat::FIMI},
});
return result;
#undef BUILD_VECTOR
}
@@ -529,10 +546,12 @@ void Benchmarks() {
<< kEol;
}
for (const auto& [dir, files, input_format] : kBenchmarks) {
GeometricMean first_cost_ratio_geomean;
GeometricMean first_speedup_geomean;
GeometricMean improvement_cost_ratio_geomean;
GeometricMean improvement_speedup_geomean;
GeometricMean element_vs_classic_cost_ratio_geomean;
GeometricMean element_vs_classic_speedup_geomean;
GeometricMean lazy_vs_classic_cost_ratio_geomean;
GeometricMean lazy_vs_classic_speedup_geomean;
GeometricMean lazy_steepest_vs_steepest_cost_ratio_geomean;
GeometricMean lazy_steepest_vs_steepest_speedup_geomean;
GeometricMean combined_cost_ratio_geomean;
GeometricMean combined_speedup_geomean;
GeometricMean tlns_cost_ratio_geomean;
@@ -570,40 +589,56 @@ void Benchmarks() {
DCHECK(classic_inv.CheckConsistency(CL::kCostAndCoverage));
}
SetCoverInvariant inv(&model);
LazyElementDegreeSolutionGenerator lazy_element_degree(&inv);
SetCoverInvariant element_degree_inv(&model);
ElementDegreeSolutionGenerator element_degree(&element_degree_inv);
CHECK(element_degree.NextSolution());
DCHECK(element_degree_inv.CheckConsistency(CL::kCostAndCoverage));
absl::StrAppend(&output, ReportBothParts(classic, element_degree));
SetCoverInvariant lazy_inv(&model);
LazyElementDegreeSolutionGenerator lazy_element_degree(&lazy_inv);
CHECK(lazy_element_degree.NextSolution());
DCHECK(inv.CheckConsistency(CL::kCostAndCoverage));
DCHECK(lazy_inv.CheckConsistency(CL::kCostAndCoverage));
absl::StrAppend(&output, ReportBothParts(classic, lazy_element_degree));
LazySteepestSearch lazy_steepest(&inv);
LazySteepestSearch lazy_steepest(&lazy_inv);
lazy_steepest.NextSolution();
absl::StrAppend(&output, ReportBothParts(steepest, lazy_steepest));
if (run_classic_solvers) {
first_cost_ratio_geomean.Add(CostRatio(classic, lazy_element_degree));
first_speedup_geomean.Add(Speedup(classic, lazy_element_degree));
improvement_cost_ratio_geomean.Add(CostRatio(steepest, lazy_steepest));
improvement_speedup_geomean.Add(Speedup(steepest, lazy_steepest));
element_vs_classic_cost_ratio_geomean.Add(
CostRatio(classic, element_degree));
element_vs_classic_speedup_geomean.Add(
Speedup(classic, element_degree));
lazy_vs_classic_cost_ratio_geomean.Add(
CostRatio(classic, lazy_element_degree));
lazy_vs_classic_speedup_geomean.Add(
Speedup(classic, lazy_element_degree));
lazy_steepest_vs_steepest_cost_ratio_geomean.Add(
CostRatio(steepest, lazy_steepest));
lazy_steepest_vs_steepest_speedup_geomean.Add(
Speedup(steepest, lazy_steepest));
combined_cost_ratio_geomean.Add(CostRatio(classic, lazy_steepest));
combined_speedup_geomean.Add(SpeedupOnCumulatedTimes(
classic, steepest, lazy_element_degree, lazy_steepest));
}
Cost best_cost = inv.cost();
SubsetBoolVector best_assignment = inv.is_selected();
Cost best_cost = lazy_inv.cost();
SubsetBoolVector best_assignment = lazy_inv.is_selected();
if (absl::GetFlag(FLAGS_tlns)) {
WallTimer timer;
LazySteepestSearch steepest(&inv);
LazySteepestSearch steepest(&lazy_inv);
timer.Start();
for (int i = 0; i < 500; ++i) {
std::vector<SubsetIndex> range =
ClearRandomSubsets(0.1 * inv.trace().size(), &inv);
ClearRandomSubsets(0.1 * lazy_inv.trace().size(), &lazy_inv);
CHECK(lazy_element_degree.NextSolution());
CHECK(steepest.NextSolution());
if (inv.cost() < best_cost) {
best_cost = inv.cost();
best_assignment = inv.is_selected();
if (lazy_inv.cost() < best_cost) {
best_cost = lazy_inv.cost();
best_assignment = lazy_inv.is_selected();
}
}
timer.Stop();
@@ -619,16 +654,21 @@ void Benchmarks() {
}
if (absl::GetFlag(FLAGS_summarize)) {
std::cout << "First speedup geometric mean: "
<< first_speedup_geomean.Get() << kEol
<< "First cost ratio geometric mean: "
<< first_cost_ratio_geomean.Get() << kEol
<< "Improvement speedup geometric mean: "
<< improvement_speedup_geomean.Get() << kEol
<< "Improvement cost ratio geometric mean: "
<< improvement_cost_ratio_geomean.Get() << kEol
<< "Combined speedup geometric mean: "
<< combined_speedup_geomean.Get() << kEol;
std::cout
<< "Element degree / classic speedup geometric mean: "
<< element_vs_classic_speedup_geomean.Get() << kEol
<< "Element degree / classic cost ratio geometric mean: "
<< element_vs_classic_cost_ratio_geomean.Get() << kEol
<< "Lazy element degree / classic speedup geometric mean: "
<< lazy_vs_classic_speedup_geomean.Get() << kEol
<< "Lazy element degree / classic cost ratio geometric mean: "
<< lazy_vs_classic_cost_ratio_geomean.Get() << kEol
<< "Improvement element degree / classic speedup geometric mean: "
<< lazy_steepest_vs_steepest_speedup_geomean.Get() << kEol
<< "Improvement cost ratio geometric mean: "
<< lazy_steepest_vs_steepest_cost_ratio_geomean.Get() << kEol
<< "Combined speedup geometric mean: "
<< combined_speedup_geomean.Get() << kEol;
if (absl::GetFlag(FLAGS_tlns)) {
std::cout << "TLNS cost ratio geometric mean: "
<< tlns_cost_ratio_geomean.Get() << kEol;

View File

@@ -0,0 +1,299 @@
// Copyright 2025 Francesco Cavaliere
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/set_cover/set_cover_submodel.h"
#include "ortools/base/stl_util.h"
namespace operations_research::scp {
template <typename SubModelT>
void PrintSubModelSummary(const SubModelT& model) {
std::cout << "SubModelView sizes: rows " << model.num_focus_elements() << "/"
<< model.num_elements() << ", columns " << model.num_focus_subsets()
<< "/" << model.num_subsets() << '\n';
std::cout << "Fixing: columns " << model.fixed_columns().size() << " cost "
<< model.fixed_cost() << '\n';
}
SubModelView::SubModelView(const Model* model)
: base_view(model, &cols_sizes_, &rows_sizes_, &cols_focus_, &rows_focus_),
full_model_(model),
cols_sizes_(model->num_subsets()),
rows_sizes_(model->num_elements()) {
DCHECK(full_model_ != nullptr);
cols_focus_.reserve(model->num_subsets());
rows_focus_.reserve(model->num_elements());
for (SubsetIndex j : model->SubsetRange()) {
cols_sizes_[j] = model->columns()[j].size();
cols_focus_.push_back(j);
}
for (ElementIndex i : model->ElementRange()) {
rows_sizes_[i] = model->rows()[i].size();
rows_focus_.push_back(i);
}
fixed_columns_.clear();
fixed_cost_ = 0.0;
PrintSubModelSummary(*this);
DCHECK(ValidateSubModel(*this));
}
SubModelView::SubModelView(const Model* model,
const std::vector<FullSubsetIndex>& columns_focus)
: base_view(model, &cols_sizes_, &rows_sizes_, &cols_focus_, &rows_focus_),
full_model_(model) {
rows_sizes_.resize(full_model_->num_elements(), 0);
for (ElementIndex i : full_model_->ElementRange()) {
rows_sizes_[i] = full_model_->rows()[i].size();
}
SetFocus(columns_focus);
}
Cost SubModelView::FixColumns(const std::vector<SubsetIndex>& columns_to_fix) {
DCHECK(full_model_ != nullptr);
Cost old_fixed_cost = fixed_cost_;
if (columns_to_fix.empty()) {
return fixed_cost_ - old_fixed_cost;
}
for (SubsetIndex j : columns_to_fix) {
DCHECK(cols_sizes_[j] > 0);
fixed_cost_ += full_model_->subset_costs()[j];
fixed_columns_.push_back(static_cast<FullSubsetIndex>(j));
cols_sizes_[j] = 0;
for (ElementIndex i : full_model_->columns()[j]) {
rows_sizes_[i] = 0;
}
}
gtl::STLEraseAllFromSequenceIf(&cols_focus_, [&](SubsetIndex j) {
DCHECK(full_model_ != nullptr);
if (cols_sizes_[j] > 0) {
cols_sizes_[j] = absl::c_count_if(full_model_->columns()[j], [&](auto i) {
return rows_sizes_[i] > 0;
});
}
return cols_sizes_[j] == 0;
});
gtl::STLEraseAllFromSequenceIf(
&rows_focus_, [&](ElementIndex i) { return !rows_sizes_[i]; });
PrintSubModelSummary(*this);
DCHECK(ValidateSubModel(*this));
return fixed_cost_ - old_fixed_cost;
}
void SubModelView::SetFocus(const std::vector<FullSubsetIndex>& columns_focus) {
DCHECK(full_model_ != nullptr);
DCHECK(!rows_sizes_.empty());
if (columns_focus.empty()) {
return;
}
cols_focus_.clear();
rows_focus_.clear();
ElementBoolVector enabled_rows(full_model_->num_elements(), false);
for (ElementIndex i : full_model_->ElementRange()) {
enabled_rows[i] = rows_sizes_[i] > 0;
}
cols_sizes_.assign(full_model_->num_subsets(), 0);
rows_sizes_.assign(full_model_->num_elements(), 0);
for (FullSubsetIndex full_j : columns_focus) {
SubsetIndex j = static_cast<SubsetIndex>(full_j);
for (ElementIndex i : full_model_->columns()[j]) {
if (enabled_rows[i] > 0) {
++cols_sizes_[j];
++rows_sizes_[i];
}
}
if (cols_sizes_[j] > 0) {
cols_focus_.push_back(j);
}
}
for (ElementIndex i : full_model_->ElementRange()) {
if (rows_sizes_[i] > 0) {
rows_focus_.push_back(i);
}
}
PrintSubModelSummary(*this);
DCHECK(ValidateSubModel(*this));
}
CoreModel::CoreModel(const Model* model)
: Model(*model),
full_model_(model),
core2full_row_map_(model->num_elements()),
full2core_row_map_(model->num_elements()),
core2full_col_map_(model->num_subsets()) {
CHECK(ElementIndex(full_model_.num_elements()) < null_element_index)
<< "Max element index is reserved.";
CHECK(SubsetIndex(full_model_.num_subsets()) < null_subset_index)
<< "Max subset index is reserved.";
absl::c_iota(core2full_row_map_, FullElementIndex());
absl::c_iota(full2core_row_map_, ElementIndex());
absl::c_iota(core2full_col_map_, FullSubsetIndex());
}
CoreModel::CoreModel(const Model* model,
const std::vector<FullSubsetIndex>& columns_focus)
: Model(),
full_model_(model),
core2full_row_map_(model->num_elements()),
full2core_row_map_(model->num_elements()) {
CHECK(ElementIndex(full_model_.num_elements()) < null_element_index)
<< "Max element index is reserved.";
CHECK(SubsetIndex(full_model_.num_subsets()) < null_subset_index)
<< "Max subset index is reserved.";
absl::c_iota(core2full_row_map_, FullElementIndex());
absl::c_iota(full2core_row_map_, ElementIndex());
SetFocus(columns_focus);
}
// Note: Assumes that columns_focus covers all rows for which rows_flags is true
// (i.e.: non-covered rows should be set to false to rows_flags). This property
// get exploited to keep the rows in the same ordering of the original model
// using "cleanish" code.
void CoreModel::SetFocus(const std::vector<FullSubsetIndex>& columns_focus) {
if (columns_focus.empty()) {
return;
}
// TODO(c4v4): change model in-place to avoid reallocations.
Model& submodel = static_cast<Model&>(*this);
submodel = Model();
core2full_col_map_.clear();
// Now we can fill the new core model
for (FullSubsetIndex full_j : columns_focus) {
core2full_col_map_.push_back(full_j);
bool first_row = true;
for (FullElementIndex full_i : full_model_.columns()[full_j]) {
ElementIndex core_i = full2core_row_map_[full_i];
if (core_i != null_element_index) {
if (first_row) {
// SetCoverModel lacks a way to remove columns
first_row = false;
submodel.AddEmptySubset(full_model_.subset_costs()[full_j]);
}
submodel.AddElementToLastSubset(core_i);
}
}
}
submodel.CreateSparseRowView();
PrintSubModelSummary(*this);
DCHECK(ValidateSubModel(*this));
}
// Mark columns and row that will be removed from the core model
// The "to-be-removed" indices are marked by setting the relative core->full
// mappings to null_*_index.
void CoreModel::MarkNewFixingInMaps(
const std::vector<SubsetIndex>& columns_to_fix) {
for (SubsetIndex old_core_j : columns_to_fix) {
fixed_cost_ += subset_costs()[old_core_j];
fixed_columns_.push_back(core2full_col_map_[old_core_j]);
core2full_col_map_[old_core_j] = null_full_subset_index;
for (ElementIndex old_core_i : columns()[old_core_j]) {
core2full_row_map_[old_core_i] = null_full_element_index;
}
}
}
// Once fixed columns and covered rows are marked, we need to create a new row
// mapping, both core->full(returned) and full->core(modifed in-place)
CoreToFullElementMapVector CoreModel::MakeOrFillBothRowMaps() {
full2core_row_map_.assign(full_model_.num_elements(), null_element_index);
CoreToFullElementMapVector new_c2f_row_map; // Future core->full mapping.
for (ElementIndex old_core_i : ElementRange()) {
FullElementIndex full_i = core2full_row_map_[old_core_i];
if (full_i != null_full_element_index) {
full2core_row_map_[full_i] = ElementIndex(new_c2f_row_map.size());
new_c2f_row_map.push_back(full_i);
}
}
return new_c2f_row_map;
}
// Create a new core model by applying the remapping from the old core model to
// the new one considering the given column fixing.
// Both the old and new core->full row mappings are required too keep track of
// what changed, the old mapping gets overwritten with the new one at the end.
// Empty columns are detected and removed - or better - not added.
Model CoreModel::MakeNewCoreModel(
const CoreToFullElementMapVector& new_c2f_row_map) {
Model new_submodel;
BaseInt new_core_j = 0;
// Loop over old core column indices.
for (SubsetIndex old_core_j : SubsetRange()) {
// If the column is not marked, then it should be mapped.
FullSubsetIndex full_j = core2full_col_map_[old_core_j];
if (full_j != null_full_subset_index) {
bool first_row = true;
// Loop over the old core column (with old core row indices).
for (ElementIndex old_core_i : columns()[old_core_j]) {
// If the row is not marked, then it should be mapped.
FullElementIndex full_i = core2full_row_map_[old_core_i];
if (full_i != null_full_element_index) {
if (first_row) {
// SetCoverModel lacks a way to remove columns
first_row = false;
new_submodel.AddEmptySubset(full_model_.subset_costs()[full_j]);
// Put the full index in the proper (new) position.
// Note that old_core_j >= new_core_j is always true.
SubsetIndex new_j(new_core_j++);
core2full_col_map_[new_j] = full_j;
}
ElementIndex new_core_i = full2core_row_map_[full_i];
DCHECK(new_core_i != null_element_index);
new_submodel.AddElementToLastSubset(new_core_i);
}
}
}
}
core2full_col_map_.resize(new_core_j);
core2full_row_map_ = std::move(new_c2f_row_map);
new_submodel.CreateSparseRowView();
return new_submodel;
}
Cost CoreModel::FixColumns(const std::vector<SubsetIndex>& columns_to_fix) {
if (columns_to_fix.empty()) {
return .0;
}
Cost old_fixed_cost = fixed_cost_;
// Mark columns to be fixed and rows that will be covered by them
MarkNewFixingInMaps(columns_to_fix);
// Compute new core->full(returned) and full->core(modified original) row maps
CoreToFullElementMapVector new_c2f_row_map = MakeOrFillBothRowMaps();
// Create new model object applying the computed mappings
static_cast<Model&>(*this) = MakeNewCoreModel(new_c2f_row_map);
PrintSubModelSummary(*this);
DCHECK(ValidateSubModel(*this));
DCHECK(absl::c_is_sorted(core2full_col_map_));
DCHECK(absl::c_is_sorted(core2full_row_map_));
return fixed_cost_ - old_fixed_cost;
}
} // namespace operations_research::scp

View File

@@ -0,0 +1,289 @@
// Copyright 2025 Francesco Cavaliere
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_SET_COVER_SET_COVER_SUBMODEL_H
#define OR_TOOLS_SET_COVER_SET_COVER_SUBMODEL_H
#include "ortools/set_cover/set_cover_views.h"
namespace operations_research::scp {
// TODO(anyone): since we are working within the scp namespace, the "SetCover*"
// prefix became redundant and can be removed. For now, only redefine it
// locally.
using Model = SetCoverModel;
// Forward declarations, see below for the definition of the classes.
struct PrimalDualState;
struct Solution;
// The CFT algorithm generates sub-models in two distinct ways:
//
// 1. It fixes specific columns (incrementally) into any generated solution.
// Once a column is fixed, it is excluded from future decisions, as it is
// already part of the solution. Additionally, rows that are covered by the
// fixed columns are removed from consideration as well, along with any
// columns that exclusively cover those rows, as they become redundant. The
// fixing process starts with the entire model and progressively fixes more
// columns until it becomes empty. A "view-based" sub-model is well-suited
// for this part.
//
// 2. The CFT mostly works on a "core" sub-model by focusing on a subset of
// columns. The core model is derived from the original model but is
// significantly smaller, as it typically includes only a limited number of
// columns per row (on average, around six columns per row). Unlike the
// incremental nature of column fixing, core models are constructed from
// scratch during each update. This type of small model can take advantage of
// a Model object which stores the sub-model explicitly in memory, avoiding
// looping over "inactive" columns and rows. Both SubModelView and CoreModel
// can be used as a core model.
//
// Two types of "core-model" representations are implemented, both of which can
// be used interchangeably:
//
// 1. SubModelView: A lightweight view of the original model. It dynamically
// filters and exposes only the active rows and columns from the original
// data structures, skipping "inactive" items.
//
// 2. CoreModel: A fully compacted and explicit representation of a sub-model.
// It stores the filtered data explicitly, making it more suitable
// for scenarios where compact storage and faster access are required.
//
// While CoreModel stores an explicit representation of the sub-model,
// SubModelView maintains vectors sized according to the original model's
// dimensions. As a result, depending on the dimensions of the original model,
// CoreModel can actually be more memory-efficient.
class SubModelView;
class CoreModel;
using SubModel = CoreModel;
// `SubModelView` provides a mechanism to interact with a subset of the rows and
// columns of a SetCoverModel, effectively creating a filtered view of the
// model. This abstraction allows operations to be performed on a restricted
// portion of the model without modifying the original data structure. The
// filtering is achieved using index lists and sizes vectors, which define the
// active rows and columns. This approach ensures flexibility and avoids
// unnecessary duplication of data. Columns/rows sizes are uses to both keep
// track of the number of elements in them and also provide the "activation"
// status: (item size == 0) <==> inactive
// SubModelView inherits from IndexListSubModelView, which provides the "view"
// machinery.
class SubModelView : public IndexListModelView {
using base_view = IndexListModelView;
public:
// Empty initialization to facilitate delayed construction
SubModelView() = default;
// Identity sub-model: all items are considered
SubModelView(const Model* model);
// Focus construction: create a sub-model with only the required items
SubModelView(const Model* model,
const std::vector<FullSubsetIndex>& columns_focus);
virtual ~SubModelView() = default;
///////// Core-model interface: /////////
// Current fixed cost: sum of the cost of the fixed columns
Cost fixed_cost() const { return fixed_cost_; }
// List of fixed columns.
const std::vector<FullSubsetIndex>& fixed_columns() const {
return fixed_columns_;
}
// Redefine the active items. The new sub-model will ignore all columns not in
// focus and (optionally) the rows for which row_flags is not true. It does
// not overwrite the current fixing.
void SetFocus(const std::vector<FullSubsetIndex>& columns_focus);
// Fix the provided columns, removing them for the submodel. Rows now covered
// by fixed columns are also removed from the submodel along with non-fixed
// columns that only cover those rows.
virtual Cost FixColumns(const std::vector<SubsetIndex>& columns_to_fix);
// Hook function for specializations. This function can be used to define a
// "small" core model considering a subset of the full model through the use
// of column-generation or by only selecting columns with good reduced cost in
// the full model.
virtual bool UpdateCore(PrimalDualState& core_state) { return false; }
private:
// Pointer to the original model
const Model* full_model_;
// Columns/rows sizes after filtering (size==0 <==> inactive)
SubsetToIntVector cols_sizes_;
ElementToIntVector rows_sizes_;
// List of columns/rows currectly active
std::vector<SubsetIndex> cols_focus_;
std::vector<ElementIndex> rows_focus_;
// Fixing data
std::vector<FullSubsetIndex> fixed_columns_;
Cost fixed_cost_ = .0;
};
// CoreModel stores a subset of the filtered columns and rows in an explicit
// Model object.
// The indices are compacted and mapped to the range [0, <sub-model-size>],
// effectively creating a smaller set-covering model. Similar to SubModelView,
// the core model supports column fixing and focusing on a subset of the
// original model. Mappings are maintained to translate indices back to the
// original model space.
class CoreModel : private Model {
public:
// Empty initialization to facilitate delayed construction
CoreModel() = default;
// Identity sub-model: all items are considered
CoreModel(const Model* model);
// Focus construction: create a sub-model with only the required items
CoreModel(const Model* model,
const std::vector<FullSubsetIndex>& columns_focus);
virtual ~CoreModel() = default;
///////// Sub-model view interface: /////////
BaseInt num_subsets() const { return full_model_.num_subsets(); }
BaseInt num_elements() const { return full_model_.num_elements(); }
BaseInt num_focus_subsets() const { return Model::num_subsets(); }
BaseInt num_focus_elements() const { return Model::num_elements(); }
BaseInt column_size(SubsetIndex j) const {
DCHECK(SubsetIndex() <= j && j < SubsetIndex(num_subsets()));
return columns()[j].size();
}
BaseInt row_size(ElementIndex i) const {
DCHECK(ElementIndex() <= i && i < ElementIndex(num_elements()));
return rows()[i].size();
}
FullElementIndex MapCoreToFullElementIndex(ElementIndex core_i) const {
DCHECK(ElementIndex() <= core_i && core_i < ElementIndex(num_elements()));
return core2full_row_map_[core_i];
}
ElementIndex MapFullToCoreElementIndex(FullElementIndex full_i) const {
DCHECK(FullElementIndex() <= full_i &&
full_i < FullElementIndex(num_elements()));
DCHECK(full2core_row_map_[full_i] != null_element_index);
return full2core_row_map_[full_i];
}
FullSubsetIndex MapCoreToFullSubsetIndex(SubsetIndex core_j) const {
DCHECK(SubsetIndex() <= core_j && core_j < SubsetIndex(num_subsets()));
return core2full_col_map_[core_j];
}
// Member function relevant for the CFT inherited from Model
using Model::columns;
using Model::ElementRange;
using Model::rows;
using Model::subset_costs;
using Model::SubsetRange;
///////// Core-model interface: /////////
// Current fixed cost: sum of the cost of the fixed columns
Cost fixed_cost() const { return fixed_cost_; }
// List of fixed columns.
const std::vector<FullSubsetIndex>& fixed_columns() const {
return fixed_columns_;
}
// Redefine the active items. The new sub-model will ignore all columns not in
// focus and (optionally) the rows for which row_flags is not true. It does
// not overwrite the current fixing.
void SetFocus(const std::vector<FullSubsetIndex>& columns_focus);
// Fix the provided columns, removing them for the submodel. Rows now covered
// by fixed columns are also removed from the submodel along with non-fixed
// columns that only cover those rows.
virtual Cost FixColumns(const std::vector<SubsetIndex>& columns_to_fix);
// Hook function for specializations. This function can be used to define a
// "small" core model considering a subset of the full model through the use
// of column-generation or by only selecting columns with good reduced cost in
// the full model.
virtual bool UpdateCore(PrimalDualState& core_state) { return false; }
private:
void MarkNewFixingInMaps(const std::vector<SubsetIndex>& columns_to_fix);
CoreToFullElementMapVector MakeOrFillBothRowMaps();
Model MakeNewCoreModel(const CoreToFullElementMapVector& new_c2f_col_map);
// Pointer to the original model
StrongModelView full_model_;
FullToCoreElementMapVector full2core_row_map_;
CoreToFullElementMapVector core2full_row_map_;
CoreToFullSubsetMapVector core2full_col_map_;
// Fixing data
Cost fixed_cost_ = .0;
std::vector<FullSubsetIndex> fixed_columns_;
static constexpr SubsetIndex null_subset_index =
std::numeric_limits<SubsetIndex>::max();
static constexpr ElementIndex null_element_index =
std::numeric_limits<ElementIndex>::max();
static constexpr FullSubsetIndex null_full_subset_index =
std::numeric_limits<FullSubsetIndex>::max();
static constexpr FullElementIndex null_full_element_index =
std::numeric_limits<FullElementIndex>::max();
};
template <typename SubModelT>
bool ValidateSubModel(const SubModelT& model) {
if (model.num_elements() <= 0) {
std::cerr << "Sub-Model has no elements.\n";
return false;
}
if (model.num_subsets() <= 0) {
std::cerr << "Sub-Model has no subsets.\n";
return false;
}
for (SubsetIndex j : model.SubsetRange()) {
const auto& column = model.columns()[j];
if (model.column_size(j) == 0) {
std::cerr << "Column " << j << " is empty.\n";
return false;
}
BaseInt j_size = std::distance(column.begin(), column.end());
if (j_size != model.column_size(j)) {
std::cerr << "Sub-Model size mismatch on column " << j << ", " << j_size
<< " != " << model.column_size(j) << "\n";
return false;
}
}
for (ElementIndex i : model.ElementRange()) {
const auto& row = model.rows()[i];
if (model.row_size(i) == 0) {
std::cerr << "Row " << i << " is empty.\n";
return false;
}
BaseInt i_size = std::distance(row.begin(), row.end());
if (i_size != model.row_size(i)) {
std::cerr << "Sub-Model size mismatch on row " << i << ", " << i_size
<< " != " << model.row_size(i) << "\n";
return false;
}
}
return true;
}
} // namespace operations_research::scp
#endif /* OR_TOOLS_SET_COVER_SET_COVER_SUBMODEL_H */

View File

@@ -488,9 +488,9 @@ TEST(SetCoverTest, KnightsCoverRandomClearMip) {
LOG(INFO) << "RandomClearMip cost: " << best_cost;
// The best solution found until 2023-08 has a cost of 350.
// http://www.contestcen.com/kn50.htm
if (BoardSize == 50) {
CHECK_GE(inv.cost(), 350);
}
// if (BoardSize == 50) {
// CHECK_GE(inv.cost(), 350);
// }
}
TEST(SetCoverTest, KnightsCoverMip) {

View File

@@ -0,0 +1,247 @@
// Copyright 2025 Francesco Cavaliere
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_SET_COVER_SET_COVER_VIEWS_H
#define OR_TOOLS_SET_COVER_SET_COVER_VIEWS_H
#include <absl/meta/type_traits.h>
#include "ortools/set_cover/base_types.h"
#include "ortools/set_cover/set_cover_model.h"
#include "views.h"
namespace operations_research {
// In the CFT algorithm, indices from different models are frequently used,
// and mixing them can lead to errors. To prevent such mistakes, strong-typed
// wrappers are employed. There are three available approaches for handling
// these indices:
// 1. Full-model strong-typed indices + {Subset,Element}Index for the core model
// 2. Core-model strong-typed indices + {Subset,Element}Index for the full model
// 3. Define new strong-typed indices both full-model and core-model
// Introducing a new set of strong-typed indices, however, can lead to a
// cascade of code duplication (or template proliferation). It also requires
// additional "view" boilerplate to properly handle the different types,
// increasing complexity.
// Currently, the simplest approach is to define only full-model indices while
// reusing the original strong types for the core model. The main challenge
// arises in FullToCoreModel, where a "filtered" full-model must be handled. In
// such cases, static casts are employed to manage the type conversions
// effectively.
DEFINE_STRONG_INT_TYPE(FullSubsetIndex, BaseInt);
DEFINE_STRONG_INT_TYPE(FullElementIndex, BaseInt);
// Syntactic sugar to define strong typed indices casts.
// Note: look at `strong_int.h` for more details about `StrongIntConvert`
#define ENABLE_EXPLICIT_STRONG_TYPE_CAST(FROM, TO) \
constexpr TO StrongIntConvert(FROM j, TO* /*unused*/) { \
return TO(static_cast<FROM::ValueType>(j)); \
}
ENABLE_EXPLICIT_STRONG_TYPE_CAST(SubsetIndex, FullSubsetIndex);
ENABLE_EXPLICIT_STRONG_TYPE_CAST(FullSubsetIndex, SubsetIndex);
ENABLE_EXPLICIT_STRONG_TYPE_CAST(ElementIndex, FullElementIndex);
ENABLE_EXPLICIT_STRONG_TYPE_CAST(FullElementIndex, ElementIndex);
using FullElementCostVector = util_intops::StrongVector<FullElementIndex, Cost>;
using FullSubsetCostVector = util_intops::StrongVector<FullSubsetIndex, Cost>;
using FullElementBoolVector = util_intops::StrongVector<FullElementIndex, bool>;
using FullSubsetBoolVector = util_intops::StrongVector<FullSubsetIndex, bool>;
// When a sub-model is created, indicies are compacted to be consecutive and
// strarting from 0 (to reduce memory usage). Core ElementIndex to original
// ElementIndex mappings are stored to translate back to the original model
// space.
using FullToCoreElementMapVector =
util_intops::StrongVector<FullElementIndex, ElementIndex>;
using CoreToFullElementMapVector =
util_intops::StrongVector<ElementIndex, FullElementIndex>;
// The same applies to SubsetIndex, which also needs to be mapped back to the
// original indexing space.
using FullToCoreSubsetMapVector =
util_intops::StrongVector<FullSubsetIndex, SubsetIndex>;
using CoreToFullSubsetMapVector =
util_intops::StrongVector<SubsetIndex, FullSubsetIndex>;
class StrongModelView {
private:
// Transformations to convert between the core and full model columns.
struct SparseColTransform {
auto operator()(const SparseColumn& column) const
-> util_intops::TransformView<
ElementIndex, ColumnEntryIndex,
util_intops::TypeCastTransform<ElementIndex, FullElementIndex>> {
return {&column};
}
};
// Transformations to convert between the core and full model rows.
struct SparseRowTransform {
auto operator()(const SparseRow& row) const -> util_intops::TransformView<
SubsetIndex, RowEntryIndex,
util_intops::TypeCastTransform<SubsetIndex, FullSubsetIndex>> {
return {&row};
}
};
public:
StrongModelView() = default;
StrongModelView(const SetCoverModel* model) : model_(model) {}
BaseInt num_subsets() const { return model_->num_subsets(); }
BaseInt num_elements() const { return model_->num_elements(); }
auto subset_costs() const
-> util_intops::TransformView<Cost, FullSubsetIndex> {
return {&model_->subset_costs()};
}
auto columns() const
-> util_intops::TransformView<SparseColumn, FullSubsetIndex,
SparseColTransform> {
return {&model_->columns()};
}
auto rows() const -> util_intops::TransformView<SparseRow, FullElementIndex,
SparseRowTransform> {
return {&model_->rows()};
}
auto SubsetRange() const -> util_intops::StrongIntRange<FullSubsetIndex> {
return {FullSubsetIndex(), FullSubsetIndex(num_subsets())};
}
auto ElementRange() const -> util_intops::StrongIntRange<FullElementIndex> {
return {FullElementIndex(), FullElementIndex(num_elements())};
}
private:
const SetCoverModel* model_;
};
class IndexListModelView {
public:
IndexListModelView() = default;
IndexListModelView(const SetCoverModel* model,
const SubsetToIntVector* cols_sizes,
const ElementToIntVector* rows_sizes,
const std::vector<SubsetIndex>* cols_focus,
const std::vector<ElementIndex>* rows_focus)
: model_(model),
cols_sizes_(cols_sizes),
rows_sizes_(rows_sizes),
cols_focus_(cols_focus),
rows_focus_(rows_focus) {}
BaseInt num_subsets() const { return model_->num_subsets(); }
BaseInt num_elements() const { return model_->num_elements(); }
BaseInt num_focus_subsets() const { return cols_focus_->size(); }
BaseInt num_focus_elements() const { return rows_focus_->size(); }
auto subset_costs() const -> util_intops::IndexListView<Cost, SubsetIndex> {
return {&model_->subset_costs(), cols_focus_};
}
auto columns() const -> util_intops::TwoLevelsView<
util_intops::IndexListView<SparseColumn, SubsetIndex>,
ElementToIntVector> {
return {{&model_->columns(), cols_focus_}, rows_sizes_};
}
auto rows() const -> util_intops::TwoLevelsView<
util_intops::IndexListView<SparseRow, ElementIndex>, SubsetToIntVector> {
return {{&model_->rows(), rows_focus_}, cols_sizes_};
}
const std::vector<SubsetIndex>& SubsetRange() const { return *cols_focus_; }
const std::vector<ElementIndex>& ElementRange() const { return *rows_focus_; }
FullElementIndex MapCoreToFullElementIndex(ElementIndex core_i) const {
DCHECK(ElementIndex() <= core_i && core_i < ElementIndex(num_elements()));
return static_cast<FullElementIndex>(core_i);
}
ElementIndex MapFullToCoreElementIndex(FullElementIndex full_i) const {
DCHECK(FullElementIndex() <= full_i &&
full_i < FullElementIndex(num_elements()));
return static_cast<ElementIndex>(full_i);
}
FullSubsetIndex MapCoreToFullSubsetIndex(SubsetIndex core_j) const {
DCHECK(SubsetIndex() <= core_j && core_j < SubsetIndex(num_subsets()));
return static_cast<FullSubsetIndex>(core_j);
}
BaseInt column_size(SubsetIndex j) const {
DCHECK(SubsetIndex() <= j && j < SubsetIndex(num_subsets()));
return (*cols_sizes_)[j];
}
BaseInt row_size(ElementIndex i) const {
DCHECK(ElementIndex() <= i && i < ElementIndex(num_elements()));
return (*rows_sizes_)[i];
}
private:
const SetCoverModel* model_;
const SubsetToIntVector* cols_sizes_;
const ElementToIntVector* rows_sizes_;
const std::vector<SubsetIndex>* cols_focus_;
const std::vector<ElementIndex>* rows_focus_;
};
// A lightweight sub-model view that uses boolean vectors to enable or disable
// specific items. Iterating over all active columns or rows is less efficient,
// particularly when only a small subset is active.
// NOTE: this view does **not** store any size-related information.
class FilterModelView {
public:
FilterModelView() = default;
FilterModelView(const SetCoverModel* model,
const SubsetBoolVector* cols_sizes,
const ElementBoolVector* rows_sizes, BaseInt num_subsets,
BaseInt num_elements)
: model_(model),
is_focus_col_(cols_sizes),
is_focus_row_(rows_sizes),
num_subsets_(num_subsets),
num_elements_(num_elements) {}
BaseInt num_subsets() const { return model_->num_subsets(); }
BaseInt num_elements() const { return model_->num_elements(); }
BaseInt num_focus_subsets() const { return num_subsets_; }
BaseInt num_focus_elements() const { return num_elements_; }
auto subset_costs() const
-> util_intops::IndexFilterView<Cost, SubsetBoolVector> {
return {&model_->subset_costs(), is_focus_col_};
}
auto columns() const -> util_intops::TwoLevelsView<
util_intops::IndexFilterView<SparseColumn, SubsetBoolVector>,
ElementBoolVector> {
return {{&model_->columns(), is_focus_col_}, is_focus_row_};
}
auto rows() const -> util_intops::TwoLevelsView<
util_intops::IndexFilterView<SparseRow, ElementBoolVector>,
SubsetBoolVector> {
return {{&model_->rows(), is_focus_row_}, is_focus_col_};
}
auto SubsetRange() const
-> util_intops::FilterIndexRangeView<SubsetIndex, SubsetBoolVector> {
return {is_focus_col_};
}
auto ElementRange() const
-> util_intops::FilterIndexRangeView<ElementIndex, ElementBoolVector> {
return {is_focus_row_};
}
private:
const SetCoverModel* model_;
const SubsetBoolVector* is_focus_col_;
const ElementBoolVector* is_focus_row_;
BaseInt num_subsets_;
BaseInt num_elements_;
};
} // namespace operations_research
#endif /* OR_TOOLS_SET_COVER_SET_COVER_VIEWS_H */

470
ortools/set_cover/views.h Normal file
View File

@@ -0,0 +1,470 @@
// Copyright 2025 Francesco Cavaliere
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_ORTOOLS_SET_COVER_VIEWS_H
#define OR_TOOLS_ORTOOLS_SET_COVER_VIEWS_H
#include <absl/meta/type_traits.h>
#include <absl/types/span.h>
#include <cstddef>
// NOTE: It may be unexpected, but views provide a subscript operator that
// directly accesses the underlying original container using the original
// indices. This behavior is particularly relevant in the context of the Set
// Cover problem, where we often work with subsets of columns or rows. Despite
// this, each column or row still contains the original indices, which are
// used for referring to other columns/rows.
//
// This design abstracts access to the underlying container, ensuring consistent
// behavior across the following scenarios:
// 1. Directly using the original container.
// 2. Accessing a subset of the original items via a view.
// 3. Using a new container with a compacted subset of items, indexing them with
// the position the have in the new container and not in the original one.
// This also need the old indices stored in rows/columns to be mapped into
// the new ones.
namespace util_intops {
namespace util {
template <typename RangeT>
using range_const_iterator_type =
absl::remove_cvref_t<decltype(std::declval<const RangeT>().begin())>;
template <typename RangeT>
using range_value_type = absl::remove_cvref_t<
decltype(*std::declval<range_const_iterator_type<RangeT>>())>;
// Enable compatibility with STL algorithms.
template <typename IterT, typename ValueT>
struct IteratorCRTP {
using iterator_category = std::input_iterator_tag;
using value_type = ValueT;
using difference_type = std::ptrdiff_t;
using pointer = IterT;
using reference = value_type&;
bool operator==(const IterT& other) const {
return !(*static_cast<const IterT*>(this) != other);
}
pointer operator->() const { return static_cast<const IterT*>(this); }
};
template <typename ValueRangeT, typename IndexT>
decltype(auto) at(const ValueRangeT* container, IndexT index) {
DCHECK(IndexT(0) <= index && index < IndexT(container->size()));
return (*container)[index];
}
} // namespace util
// View exposing only the indices of a container that are marked as active.
// Looping over this view is equivalent to:
//
// for (IntT index; index < IntT<is_active_->size()); ++index) {
// if (is_active[index]) {
// your_code(index);
// }
// }
//
template <typename IntT, typename EnableVectorT>
class FilterIndexRangeView {
public:
struct FilterIndicesViewIterator
: util::IteratorCRTP<FilterIndicesViewIterator, IntT> {
FilterIndicesViewIterator(IntT index, const EnableVectorT* is_active)
: index_(index), is_active_(is_active) {
AdjustToValidValue();
std::vector<bool> vb;
}
bool operator!=(FilterIndicesViewIterator other) const {
return index_ != other.index_;
}
FilterIndicesViewIterator& operator++() {
++index_;
AdjustToValidValue();
return *this;
}
IntT operator*() const { return index_; }
private:
void AdjustToValidValue() {
while (index_ < IntT(is_active_->size()) &&
!util::at(is_active_, index_)) {
++index_;
}
}
IntT index_;
const EnableVectorT* is_active_;
};
FilterIndexRangeView(const EnableVectorT* is_active)
: is_active_(is_active) {}
FilterIndicesViewIterator begin() const {
return FilterIndicesViewIterator(IntT(0), is_active_);
}
FilterIndicesViewIterator end() const {
return FilterIndicesViewIterator(IntT(is_active_->size()), is_active_);
}
private:
const EnableVectorT* is_active_;
};
// View exposing only the elements of a container that are indexed by a list of
// indices.
// Looping over this view is equivalent to:
//
// for (decltype(auto) index : indices) {
// your_code(container[index]);
// }
//
template <typename ValueT, typename IndexT>
class IndexListView {
public:
using value_type = const ValueT;
using index_type = const IndexT;
using value_iterator = typename absl::Span<value_type>::iterator;
using index_iterator = typename absl::Span<index_type>::iterator;
struct IndexListViewIterator
: util::IteratorCRTP<IndexListViewIterator, value_type> {
IndexListViewIterator(absl::Span<value_type> values, index_iterator iter)
: values_(values), index_iter_(iter) {}
bool operator!=(const IndexListViewIterator& other) const {
DCHECK_EQ(values_.data(), other.values_.data());
return index_iter_ != other.index_iter_;
}
IndexListViewIterator& operator++() {
++index_iter_;
return *this;
}
decltype(auto) operator*() const {
return values_[static_cast<size_t>(*index_iter_)];
}
private:
absl::Span<value_type> values_;
index_iterator index_iter_;
};
IndexListView() = default;
template <typename ValueRangeT, typename IndexRangeT>
IndexListView(const ValueRangeT* values, const IndexRangeT* indices)
: values_(absl::MakeConstSpan(*values)),
indices_(absl::MakeConstSpan(*indices)) {}
auto size() const { return indices_.size(); }
bool empty() const { return indices_.empty(); }
// NOTE: uses indices of the original container, not the filtered one
decltype(auto) operator[](index_type index) const {
// TODO(user): we could check that index is in indices_, but that's is O(n).
return values_[static_cast<size_t>(index)];
}
IndexListViewIterator begin() const {
return IndexListViewIterator(values_, indices_.begin());
}
IndexListViewIterator end() const {
return IndexListViewIterator(values_, indices_.end());
}
private:
absl::Span<value_type> values_;
absl::Span<index_type> indices_;
};
// This view provides access to elements of a container of integral types, which
// are used as indices to a filter (in) vector.
// The view filters and returns only those elements that, when used a subscript
// to the filter vector, are evaluated to a true value.
// Looping over this view is equivalent to:
//
// for (decltype(auto) item : container) {
// if (is_active[item]) {
// your_code(iter);
// }
// }
//
template <typename ValueT, typename EnableVectorT>
class ValueFilterView {
public:
using value_type = const ValueT;
using value_iterator = typename absl::Span<value_type>::iterator;
struct ValueFilterViewIterator
: util::IteratorCRTP<ValueFilterViewIterator, value_type> {
ValueFilterViewIterator(value_iterator iterator, value_iterator end,
const EnableVectorT* is_active)
: iterator_(iterator), end_(end), is_active_(is_active) {
AdjustToValidValue();
}
bool operator!=(const ValueFilterViewIterator& other) const {
DCHECK_EQ(is_active_, other.is_active_);
return iterator_ != other.iterator_;
}
ValueFilterViewIterator& operator++() {
++iterator_;
AdjustToValidValue();
return *this;
}
decltype(auto) operator*() const { return *iterator_; }
private:
void AdjustToValidValue() {
while (iterator_ != end_ && !util::at(is_active_, *iterator_)) {
++iterator_;
}
}
value_iterator iterator_;
value_iterator end_;
const EnableVectorT* is_active_;
};
template <typename ValueRangeT>
ValueFilterView(const ValueRangeT* values, const EnableVectorT* is_active)
: values_(absl::MakeConstSpan(*values)), is_active_(is_active) {
DCHECK(values != nullptr);
DCHECK(is_active != nullptr);
}
ValueFilterViewIterator begin() const {
return ValueFilterViewIterator(values_.begin(), values_.end(), is_active_);
}
ValueFilterViewIterator end() const {
return ValueFilterViewIterator(values_.end(), values_.end(), is_active_);
}
// NOTE: uses indices of the original container, not the filtered one
template <typename IndexT>
decltype(auto) operator[](IndexT index) const {
decltype(auto) value = values_[static_cast<size_t>(index)];
DCHECK(util::at(is_active_, value))
<< "Inactive value. Are you using relative indices?";
return value;
}
private:
absl::Span<value_type> values_;
const EnableVectorT* is_active_;
};
// Somewhat equivalent to ValueFilterView<StrongIntRange, EnableVectorT>
// Looping over this view is equivalent to:
//
// auto c_it = container.begin();
// auto active_it = is_active.begin();
// for (; active_it != is_active.end(); ++active_it, ++c_it) {
// if (*active_it) {
// your_code(*c_it);
// }
// }
//
template <typename ValueT, typename EnableVectorT>
class IndexFilterView {
public:
using value_type = const ValueT;
using value_iterator = typename absl::Span<value_type>::iterator;
using enable_iterator = util::range_const_iterator_type<EnableVectorT>;
struct IndexFilterViewIterator
: util::IteratorCRTP<IndexFilterViewIterator, value_type> {
IndexFilterViewIterator(value_iterator iterator,
enable_iterator is_active_begin,
enable_iterator is_active_end)
: iterator_(iterator),
is_active_iter_(is_active_begin),
is_active_end_(is_active_end) {
AdjustToValidValue();
}
bool operator!=(const IndexFilterViewIterator& other) const {
DCHECK(is_active_end_ == other.is_active_end_);
return iterator_ != other.iterator_;
}
IndexFilterViewIterator& operator++() {
++is_active_iter_;
++iterator_;
AdjustToValidValue();
return *this;
}
decltype(auto) operator*() const { return *iterator_; }
private:
void AdjustToValidValue() {
while (is_active_iter_ != is_active_end_ && !*is_active_iter_) {
++is_active_iter_;
++iterator_;
}
}
value_iterator iterator_;
enable_iterator is_active_iter_;
enable_iterator is_active_end_;
};
template <typename ValueRangeT>
IndexFilterView(const ValueRangeT* values, const EnableVectorT* is_active_)
: values_(absl::MakeConstSpan(*values)), is_active_(is_active_) {
DCHECK(values != nullptr);
DCHECK(is_active_ != nullptr);
DCHECK_EQ(values->size(), is_active_->size());
}
IndexFilterViewIterator begin() const {
return IndexFilterViewIterator(values_.begin(), is_active_->begin(),
is_active_->end());
}
IndexFilterViewIterator end() const {
return IndexFilterViewIterator(values_.end(), is_active_->end(),
is_active_->end());
}
// NOTE: uses indices of the original container, not the filtered one
template <typename IndexT>
decltype(auto) operator[](IndexT index) const {
DCHECK(util::at(is_active_, index))
<< "Inactive value. Are you using relative indices?";
return values_[static_cast<size_t>(index)];
}
private:
absl::Span<value_type> values_;
const EnableVectorT* is_active_;
};
// This view provides a mechanism to access and filter elements in a 2D
// container. The filtering is applied in two stages:
// 1. The first dimension is generic and can be both and index-list or
// bool-vector based view.
// 2. The second dimension (items of each sub-container) is further filtered
// using a boolean-vector-like object, which determines which elements within
// the sub-container are included in the view.
template <typename Lvl1ViewT, typename EnableVectorT>
class TwoLevelsView : public Lvl1ViewT {
public:
using level1_iterator = util::range_const_iterator_type<Lvl1ViewT>;
using level1_value = util::range_value_type<Lvl1ViewT>;
using level2_value = util::range_value_type<level1_value>;
using level2_type = ValueFilterView<level2_value, EnableVectorT>;
struct TwoLevelsViewIterator
: util::IteratorCRTP<TwoLevelsViewIterator, level2_type> {
TwoLevelsViewIterator(level1_iterator iter,
const EnableVectorT* active_items)
: iter_(iter), active_items_(active_items) {}
bool operator!=(const TwoLevelsViewIterator& other) const {
DCHECK_EQ(active_items_, other.active_items_);
return iter_ != other.iter_;
}
TwoLevelsViewIterator& operator++() {
++iter_;
return *this;
}
level2_type operator*() const {
return level2_type(&(*iter_), active_items_);
}
private:
level1_iterator iter_;
const EnableVectorT* active_items_;
};
TwoLevelsView() = default;
TwoLevelsView(Lvl1ViewT lvl1_view, const EnableVectorT* active_items)
: Lvl1ViewT(lvl1_view), active_items_(active_items) {}
TwoLevelsViewIterator begin() const {
return TwoLevelsViewIterator(Lvl1ViewT::begin(), active_items_);
}
TwoLevelsViewIterator end() const {
return TwoLevelsViewIterator(Lvl1ViewT::end(), active_items_);
}
template <typename indexT>
level2_type operator[](indexT i) const {
auto& level2_container = Lvl1ViewT::operator[](i);
return level2_type(&level2_container, active_items_);
}
private:
const EnableVectorT* active_items_;
};
struct NoTransform {
template <typename T>
T&& operator()(T&& value) const {
return std::forward<T>(value);
}
};
template <typename FromT, typename ToT>
struct TypeCastTransform {
ToT operator()(FromT value) const { return static_cast<ToT>(value); }
};
// View applying stateless transformation to the values of a continugous
// container. Looping over this view is equivalent to:
//
// for (IndexT i(0); i < IndexT(container.size()); ++i) {
// your_code(ValueTransformT()(values_[static_cast<size_t>(i)]));
// }
//
template <typename ValueT, typename IndexT,
typename ValueTransformT = NoTransform>
class TransformView {
public:
using value_type = const ValueT;
using value_iterator = typename absl::Span<value_type>::iterator;
struct TransformViewIterator
: util::IteratorCRTP<TransformViewIterator, value_type> {
TransformViewIterator(value_iterator iterator) : iterator_(iterator) {}
bool operator!=(const TransformViewIterator& other) const {
return iterator_ != other.iterator_;
}
TransformViewIterator& operator++() {
++iterator_;
return *this;
}
decltype(auto) operator*() const { return ValueTransformT()(*iterator_); }
private:
value_iterator iterator_;
};
TransformView() = default;
template <typename ValueRangeT>
TransformView(const ValueRangeT* values)
: values_(absl::MakeConstSpan(*values)) {}
auto size() const { return values_.size(); }
bool empty() const { return values_.empty(); }
decltype(auto) operator[](IndexT index) const {
return ValueTransformT()(values_[static_cast<size_t>(index)]);
}
TransformViewIterator begin() const {
return TransformViewIterator(values_.begin());
}
TransformViewIterator end() const {
return TransformViewIterator(values_.end());
}
private:
absl::Span<value_type> values_;
};
} // namespace util_intops
#endif /* OR_TOOLS_ORTOOLS_SET_COVER_VIEWS_H */