always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'main' of github.com:google/or-tools

Browse Source
This commit is contained in:
Laurent Perron
2023-12-05 22:17:10 +01:00
parent 7ee0f88644 43bbadb525
commit 308b43552b
24 changed files with 610 additions and 152 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
ortools/algorithms/BUILD.bazel
View File

@@ -99,7 +99,6 @@ cc_library(
)
# Weighted set covering
cc_library(
name = "set_cover_model",
srcs = ["set_cover_model.cc"],
@@ -123,6 +122,18 @@ cc_library(
],
)
cc_library(
name = "set_cover_mip",
srcs = ["set_cover_mip.cc"],
hdrs = ["set_cover_mip.h"],
deps = [
":set_cover_ledger",
"//ortools/linear_solver",
"//ortools/lp_data:base",
"@com_google_absl//absl/log:check",
],
)
cc_library(
name = "set_cover_utils",
srcs = ["set_cover_utils.cc"],
@@ -139,6 +150,7 @@ cc_library(
srcs = ["set_cover.cc"],
hdrs = ["set_cover.h"],
deps = [
":set_cover_mip",
":set_cover_utils",
"//ortools/base",
"@com_google_absl//absl/log",

70
ortools/algorithms/README.md Normal file
View File

@@ -0,0 +1,70 @@
# Various algorithms
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.h`](set_cover.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();
```

95
ortools/algorithms/set_cover.cc
View File

@@ -14,10 +14,13 @@
#include "ortools/algorithms/set_cover.h"
#include <algorithm>
#include <cstddef>
#include <iterator>
#include <limits>
#include <numeric>
#include <vector>
#include "absl/container/flat_hash_set.h"
#include "absl/log/check.h"
#include "absl/random/random.h"
#include "ortools/algorithms/set_cover_ledger.h"
@@ -289,31 +292,97 @@ bool GuidedTabuSearch::NextSolution(const std::vector<SubsetIndex>& focus,
return true;
}
std::vector<SubsetIndex> ClearProportionRandomly(double proportion,
namespace {
void SampleSubsets(std::vector<SubsetIndex>* list, std::size_t num_subsets) {
num_subsets = std::min(num_subsets, list->size());
CHECK_GE(num_subsets, 0);
std::shuffle(list->begin(), list->end(), absl::BitGen());
list->resize(num_subsets);
}
} // namespace
std::vector<SubsetIndex> ClearRandomSubsets(std::size_t num_subsets,
SetCoverLedger* ledger) {
return ClearProportionRandomly(ledger->model()->all_subsets(), proportion,
return ClearRandomSubsets(ledger->model()->all_subsets(), num_subsets,
ledger);
}
std::vector<SubsetIndex> ClearProportionRandomly(
const std::vector<SubsetIndex>& focus, double proportion,
std::vector<SubsetIndex> ClearRandomSubsets(
const std::vector<SubsetIndex>& focus, std::size_t num_subsets,
SetCoverLedger* ledger) {
CHECK_LE(proportion, 1.0);
CHECK_GE(proportion, 0.0);
std::vector<SubsetIndex> choice_indices;
num_subsets = std::min(num_subsets, focus.size());
CHECK_GE(num_subsets, 0);
std::vector<SubsetIndex> chosen_indices;
for (const SubsetIndex subset : focus) {
if (ledger->is_selected(subset)) {
choice_indices.push_back(subset);
chosen_indices.push_back(subset);
}
}
std::shuffle(choice_indices.begin(), choice_indices.end(), absl::BitGen());
const int num_subsets_to_clear = proportion * choice_indices.size();
choice_indices.resize(num_subsets_to_clear);
for (const SubsetIndex subset : choice_indices) {
SampleSubsets(&chosen_indices, num_subsets);
for (const SubsetIndex subset : chosen_indices) {
// Use UnsafeToggle because we allow non-solutions.
ledger->UnsafeToggle(subset, false);
}
return choice_indices;
return chosen_indices;
}
std::vector<SubsetIndex> ClearMostCoveredElements(std::size_t num_subsets,
SetCoverLedger* ledger) {
return ClearMostCoveredElements(ledger->model()->all_subsets(), num_subsets,
ledger);
}
std::vector<SubsetIndex> ClearMostCoveredElements(
const std::vector<SubsetIndex>& focus, std::size_t num_subsets,
SetCoverLedger* ledger) {
// This is the vector we will return.
std::vector<SubsetIndex> chosen_indices;
const ElementToSubsetVector& coverage = ledger->coverage();
// Compute a permutation of the element indices by decreasing order of
// coverage by element.
std::vector<ElementIndex> permutation(coverage.size().value());
std::iota(permutation.begin(), permutation.end(), 0);
std::sort(permutation.begin(), permutation.end(),
[&coverage](ElementIndex i, ElementIndex j) {
return coverage[i] > coverage[j];
});
// Now, for the elements that are over-covered (coverage > 1), collect the
// sets that are used.
absl::flat_hash_set<SubsetIndex> used_subsets_collection;
for (ElementIndex element : permutation) {
if (coverage[element] <= 1) break;
for (SubsetIndex subset : ledger->model()->rows()[element]) {
if (ledger->is_selected(subset)) {
used_subsets_collection.insert(subset);
}
}
}
// Now the impacted subset is a vector representation of the flat_hash_set
// collection.
std::vector<SubsetIndex> impacted_subsets(used_subsets_collection.begin(),
used_subsets_collection.end());
// Sort the impacted subsets to be able to intersect the vector later.
std::sort(impacted_subsets.begin(), impacted_subsets.end());
// chosen_indices = focus ⋂ impacted_subsets
std::set_intersection(focus.begin(), focus.end(), impacted_subsets.begin(),
impacted_subsets.end(),
std::back_inserter(chosen_indices));
std::shuffle(chosen_indices.begin(), chosen_indices.end(), absl::BitGen());
chosen_indices.resize(std::min(chosen_indices.size(), num_subsets));
// Sort before traversing indices (and memory) in order.
std::sort(chosen_indices.begin(), chosen_indices.end());
for (const SubsetIndex subset : chosen_indices) {
// Use UnsafeToggle because we allow non-solutions.
ledger->UnsafeToggle(subset, false);
}
return chosen_indices;
}
} // namespace operations_research

34
ortools/algorithms/set_cover.h
View File

@@ -14,9 +14,11 @@
#ifndef OR_TOOLS_ALGORITHMS_SET_COVER_H_
#define OR_TOOLS_ALGORITHMS_SET_COVER_H_
#include <cstddef>
#include <vector>
#include "ortools/algorithms/set_cover_ledger.h"
#include "ortools/algorithms/set_cover_model.h"
#include "ortools/algorithms/set_cover_utils.h"
namespace operations_research {
@@ -94,13 +96,19 @@ class RandomSolutionGenerator {
// remaining uncovered elements as possible for the least possible cost per
// element and iterate.
//
// The following paper dives into the details of this class of algorithms.
// The following papers dive into the details of this class of algorithms.
//
// Neal E. Young, Greedy Set-Cover Algorithms (1974-1979, Chvatal,
// Johnson, Lovasz, Stein), in Kao, ed., Encyclopedia of Algorithms. Draft at:
// http://www.cs.ucr.edu/~neal/non_arxiv/Young08SetCover.pdf
//
// G. Cormode, H. Karloff, A. Wirth (2010) "Set Cover Algorithms for Very Large
// Datasets", CIKM'10, ACM, 2010.
//
// Tommaso Adamo, Gianpaolo Ghiani, Emanuela Guerriero and Deborah Pareo "A
// Surprisal-Based Greedy Heuristic for the Set Covering Problem". Algorithms,
// 2023.
class GreedySolutionGenerator {
public:
explicit GreedySolutionGenerator(SetCoverLedger* ledger)
@@ -258,13 +266,27 @@ class GuidedTabuSearch {
TabuList<SubsetIndex> tabu_list_;
};
// Randomly clears a proportion (between 0 and 1) of the solution.
std::vector<SubsetIndex> ClearProportionRandomly(double proportion,
SetCoverLedger* ledger);
// Randomly clears a proportion num_subsets variables in the solution.
// Returns a list of subset indices to be potentially reused as a focus.
std::vector<SubsetIndex> ClearRandomSubsets(std::size_t num_subsets,
SetCoverLedger* ledger);
std::vector<SubsetIndex> ClearProportionRandomly(
const std::vector<SubsetIndex>& focus, double proportion,
// Same as above, but clears the subset indices in focus.
std::vector<SubsetIndex> ClearRandomSubsets(
const std::vector<SubsetIndex>& focus, std::size_t num_subsets,
SetCoverLedger* ledger);
// Clears the variables that cover the most covered elements. This is capped
// by num_subsets.
// Return the list of chosen subset indices to be potentially reused as a focus.
std::vector<SubsetIndex> ClearMostCoveredElements(std::size_t num_subsets,
SetCoverLedger* ledger);
// Same as above, but clears the subset indices in focus.
std::vector<SubsetIndex> ClearMostCoveredElements(
const std::vector<SubsetIndex>& focus, std::size_t num_subsets,
SetCoverLedger* ledger);
} // namespace operations_research
#endif // OR_TOOLS_ALGORITHMS_SET_COVER_H_

106
ortools/algorithms/set_cover_mip.cc Normal file
View File

@@ -0,0 +1,106 @@
// Copyright 2010-2022 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/algorithms/set_cover_mip.h"
#include <cstdint>
#include <limits>
#include <vector>
#include "ortools/algorithms/set_cover_ledger.h"
#include "ortools/algorithms/set_cover_model.h"
#include "ortools/base/logging.h"
#include "ortools/linear_solver/linear_solver.h"
#include "ortools/lp_data/lp_types.h"
namespace operations_research {
template <typename IndexType, typename ValueType>
using StrictVector = glop::StrictITIVector<IndexType, ValueType>;
bool SetCoverMip::NextSolution() {
return NextSolution(ledger_->model()->all_subsets());
}
bool SetCoverMip::NextSolution(const std::vector<SubsetIndex>& focus) {
SetCoverModel* model = ledger_->model();
const SubsetIndex num_subsets(model->num_subsets());
const ElementIndex num_elements(model->num_elements());
SubsetBoolVector choices = ledger_->GetSolution();
MPSolver::OptimizationProblemType problem_type;
switch (mip_solver_) {
case SetCoverMipSolver::SCIP:
problem_type = MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING;
break;
case SetCoverMipSolver::GUROBI:
problem_type = MPSolver::GUROBI_MIXED_INTEGER_PROGRAMMING;
break;
case SetCoverMipSolver::SAT:
problem_type = MPSolver::SAT_INTEGER_PROGRAMMING;
break;
default:
LOG(WARNING) << "Unknown solver value, defaulting to SCIP";
problem_type = MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING;
}
// We are using MPSolver, which is deprecated, because MathOpt does not
// provide an interface without using protobufs.
// We construct a restricted MIP, omitting all the parts of the problem
// that are already fixed in the ledger. The goal is to not spend time sending
// data, and having the MIP solver re-discover the fixed variables.
MPSolver solver("set cover mip", problem_type);
solver.SuppressOutput();
MPObjective* const objective = solver.MutableObjective();
objective->SetMinimization();
StrictVector<ElementIndex, MPConstraint*> constraints(num_elements, nullptr);
StrictVector<SubsetIndex, MPVariable*> vars(num_subsets, nullptr);
for (const SubsetIndex subset : focus) {
vars[subset] = solver.MakeBoolVar("");
objective->SetCoefficient(vars[subset], model->subset_costs()[subset]);
for (ElementIndex element : model->columns()[subset]) {
if (ledger_->coverage(element) > 0) continue;
if (constraints[element] == nullptr) {
constexpr double kInfinity = std::numeric_limits<double>::infinity();
constraints[element] = solver.MakeRowConstraint(1.0, kInfinity);
}
constraints[element]->SetCoefficient(vars[subset], 1);
}
}
// set_time_limit takes milliseconds as a unit.
solver.set_time_limit(static_cast<int64_t>(time_limit_in_seconds_ * 1000));
// Call the solver.
const MPSolver::ResultStatus solve_status = solver.Solve();
switch (solve_status) {
case MPSolver::OPTIMAL:
break;
case MPSolver::FEASIBLE:
break;
case MPSolver::INFEASIBLE:
LOG(ERROR) << "Did not find solution. Problem is infeasible.";
break;
case MPSolver::UNBOUNDED:
LOG(ERROR) << "Did not find solution. Problem is unbounded.";
break;
default:
LOG(ERROR) << "Solving resulted in an error.";
return false;
}
for (const SubsetIndex subset : focus) {
choices[subset] = (vars[subset]->solution_value() > 0.9);
}
ledger_->LoadSolution(choices);
return true;
}
} // namespace operations_research

58
ortools/algorithms/set_cover_mip.h Normal file
View File

@@ -0,0 +1,58 @@
// Copyright 2010-2022 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.
#ifndef OR_TOOLS_ALGORITHMS_SET_COVER_MIP_H_
#define OR_TOOLS_ALGORITHMS_SET_COVER_MIP_H_
#include <vector>
#include "ortools/algorithms/set_cover_ledger.h"
#include "ortools/algorithms/set_cover_model.h"
namespace operations_research {
enum class SetCoverMipSolver : int { SCIP = 0, SAT = 1, GUROBI = 2 };
class SetCoverMip {
public:
explicit SetCoverMip(SetCoverLedger* ledger)
: ledger_(ledger),
mip_solver_(SetCoverMipSolver::SCIP),
time_limit_in_seconds_(0.02) {}
// Returns true if a solution was found.
// TODO(user): Add time-outs and exit with a partial solution. This seems
// unlikely, though.
bool NextSolution();
// Computes the next partial solution considering only the subsets whose
// indices are in focus.
bool NextSolution(const std::vector<SubsetIndex>& focus);
void SetMipSolver(const SetCoverMipSolver mip_solver) {
mip_solver_ = mip_solver;
}
void SetTimeLimitInSeconds(double limit) { time_limit_in_seconds_ = limit; }
private:
// The ledger on which the algorithm will run.
SetCoverLedger* ledger_;
// The MIP solver flavor used by the instance.
SetCoverMipSolver mip_solver_;
double time_limit_in_seconds_;
};
} // namespace operations_research
#endif // OR_TOOLS_ALGORITHMS_SET_COVER_MIP_H_

96
ortools/algorithms/set_cover_test.cc
View File

@@ -15,13 +15,13 @@
#include <limits>
#include <string>
#include <vector>
#include "absl/log/check.h"
#include "absl/strings/str_cat.h"
#include "benchmark/benchmark.h"
#include "gtest/gtest.h"
#include "ortools/algorithms/set_cover_ledger.h"
#include "ortools/algorithms/set_cover_mip.h"
#include "ortools/algorithms/set_cover_model.h"
#include "ortools/base/logging.h"
@@ -187,7 +187,11 @@ TEST(SetCoverTest, KnightsCoverTrivial) {
}
TEST(SetCoverTest, KnightsCoverGreedyAndTabu) {
const int BoardSize = 15;
#ifdef NDEBUG
constexpr int BoardSize = 50;
#else
constexpr int BoardSize = 15;
#endif
SetCoverModel model = CreateKnightsCoverModel(BoardSize, BoardSize);
SetCoverLedger ledger(&model);
@@ -208,37 +212,99 @@ TEST(SetCoverTest, KnightsCoverGreedyAndTabu) {
}
TEST(SetCoverTest, KnightsCoverRandomClear) {
const int BoardSize = 15;
#ifdef NDEBUG
constexpr int BoardSize = 50;
#else
constexpr int BoardSize = 15;
#endif
SetCoverModel model = CreateKnightsCoverModel(BoardSize, BoardSize);
SetCoverLedger ledger(&model);
Cost best_cost = std::numeric_limits<Cost>::max();
SubsetBoolVector best_choices = ledger.GetSolution();
std::vector<SubsetIndex> focus = model.all_subsets();
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(focus));
// LOG(INFO) << "GreedySolutionGenerator cost: " << ledger.cost();
CHECK(greedy.NextSolution());
SteepestSearch steepest(&ledger);
CHECK(steepest.NextSolution(focus, 10000));
// LOG(INFO) << "SteepestSearch cost: " << ledger.cost();
CHECK(steepest.NextSolution(10000));
EXPECT_TRUE(ledger.CheckSolution());
if (ledger.cost() < best_cost) {
best_cost = ledger.cost();
best_choices = ledger.GetSolution();
LOG(INFO) << "Best cost: " << best_cost << " at iteration = " << i;
}
ledger.LoadSolution(best_choices);
ClearProportionRandomly(0.1, &ledger);
// focus = ledger.ComputeSettableSubsets();
}
if (ledger.cost() < 0) {
// The best solution found until now has a cost of 350.
// http://www.contestcen.com/kn50.htm
DisplayKnightsCoverSolution(best_choices, BoardSize, BoardSize);
ledger.LoadSolution(best_choices);
DisplayKnightsCoverSolution(best_choices, BoardSize, BoardSize);
LOG(INFO) << "RandomClear 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(ledger.cost(), 350);
}
}
TEST(SetCoverTest, KnightsCoverRandomClearMip) {
#ifdef NDEBUG
constexpr int BoardSize = 50;
#else
constexpr int BoardSize = 15;
#endif
SetCoverModel model = CreateKnightsCoverModel(BoardSize, BoardSize);
SetCoverLedger ledger(&model);
Cost best_cost = std::numeric_limits<Cost>::max();
SubsetBoolVector best_choices = ledger.GetSolution();
GreedySolutionGenerator greedy(&ledger);
CHECK(greedy.NextSolution());
LOG(INFO) << "GreedySolutionGenerator cost: " << ledger.cost();
SteepestSearch steepest(&ledger);
CHECK(steepest.NextSolution(10000));
LOG(INFO) << "SteepestSearch cost: " << ledger.cost();
best_cost = ledger.cost();
best_choices = ledger.GetSolution();
for (int i = 0; i < 100; ++i) {
ledger.LoadSolution(best_choices);
auto focus = ClearRandomSubsets(0.1 * model.num_subsets().value(), &ledger);
SetCoverMip mip(&ledger);
mip.SetTimeLimitInSeconds(1);
mip.NextSolution(focus);
EXPECT_TRUE(ledger.CheckSolution());
if (ledger.cost() < best_cost) {
best_cost = ledger.cost();
best_choices = ledger.GetSolution();
LOG(INFO) << "Best cost: " << best_cost << " at iteration = " << i;
}
}
ledger.LoadSolution(best_choices);
DisplayKnightsCoverSolution(best_choices, BoardSize, BoardSize);
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(ledger.cost(), 350);
}
}
TEST(SetCoverTest, KnightsCoverMip) {
#ifdef NDEBUG
constexpr int BoardSize = 50;
#else
constexpr int BoardSize = 15;
#endif
SetCoverModel model = CreateKnightsCoverModel(BoardSize, BoardSize);
SetCoverLedger ledger(&model);
SetCoverMip mip(&ledger);
mip.SetTimeLimitInSeconds(10);
mip.NextSolution();
SubsetBoolVector best_choices = ledger.GetSolution();
DisplayKnightsCoverSolution(best_choices, BoardSize, BoardSize);
LOG(INFO) << "Mip cost: " << ledger.cost();
}
void BM_Steepest(benchmark::State& state) {

9
ortools/glop/lu_factorization.cc
View File

@@ -18,6 +18,7 @@
#include <cstdlib>
#include <vector>
#include "absl/log/check.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/lp_data/lp_utils.h"
@@ -65,9 +66,11 @@ Status LuFactorization::ComputeFactorization(
stats_.basis_num_entries.Add(matrix.num_entries().value());
});
// TODO(user): This might fail on badly scaled matrices.
// I still prefer to keep it as a DCHECK for tests though.
DCHECK(CheckFactorization(matrix, Fractional(1e-6)));
// Note(user): This might fail on badly scaled matrices. I still prefer to
// keep it as a DCHECK() for tests though, but I only test it for small
// matrices.
DCHECK(matrix.num_rows() > 100 ||
CheckFactorization(matrix, Fractional(1e-6)));
return Status::OK();
}

6
ortools/glop/pricing.h
View File

@@ -14,9 +14,11 @@
#ifndef OR_TOOLS_GLOP_PRICING_H_
#define OR_TOOLS_GLOP_PRICING_H_
#include <cmath>
#include <random>
#include <string>
#include "absl/log/check.h"
#include "absl/random/bit_gen_ref.h"
#include "absl/random/random.h"
#include "ortools/lp_data/lp_types.h"
@@ -178,7 +180,7 @@ inline void DynamicMaximum<Index>::StartDenseUpdates() {
template <typename Index>
inline void DynamicMaximum<Index>::DenseAddOrUpdate(Index position,
Fractional value) {
DCHECK(IsFinite(value));
DCHECK(!std::isnan(value));
DCHECK(tops_.empty());
is_candidate_.Set(position);
values_[position] = value;
@@ -187,7 +189,7 @@ inline void DynamicMaximum<Index>::DenseAddOrUpdate(Index position,
template <typename Index>
inline void DynamicMaximum<Index>::AddOrUpdate(Index position,
Fractional value) {
DCHECK(IsFinite(value));
DCHECK(!std::isnan(value));
is_candidate_.Set(position);
values_[position] = value;
if (value >= threshold_) UpdateTopK(position, value);

9
ortools/graph/README.md
View File

@@ -4,9 +4,6 @@ This directory contains data structures and algorithms for graph and
network flow problems.
It contains in particular:
* a compact and efficient graph representation,
[`EbertGraph`](https://dl.acm.org/doi/abs/10.1145/214762.214769),
which is used for most of the algorithms herein, unless specified otherwise.
* well-tuned algorithms (for example shortest, paths and
[Hamiltonian paths](https://en.wikipedia.org/wiki/Hamiltonian_path)).
* hard-to-find algorithms (Hamiltonian paths, push-relabel flow algorithms).
@@ -14,8 +11,7 @@ It contains in particular:
Graph representations:
* [ebert_graph.h](./ebert_graph.h): entry point for a directed graph class.
* [digraph.h](./digraph.h): entry point for a directed graph class.
To be deprecated by `ebert_graph.h`.
Deprecated. Prefer using [`//util/graph/graph.h`](../../graph/graph.h).
Paths:
@@ -37,8 +33,7 @@ Graph decompositions:
or `digraph.h`.)
* [strongly_connected_components.h](./strongly_connected_components.h): entry
point for computing the strongly connected components of a directed graph,
based on Tarjan's algorithm.
point for computing the strongly connected components of a directed graph.
* [cliques.h](./cliques.h): entry point for computing maximum cliques and
clique covers in a directed graph, based on the Bron-Kerbosch algorithm.

30
ortools/linear_solver/gurobi_interface.cc
View File

@@ -236,7 +236,11 @@ class GurobiInterface : public MPSolverInterface {
int SolutionCount() const;
GRBmodel* model_;
GRBenv* env_;
// The environment used to create model_. Note that once the model is created,
// it used a copy of the environment, accessible via GRBgetenv(model_), which
// you should use for setting parameters. Use global_env_ only to create a new
// model or for GRBgeterrormsg().
GRBenv* global_env_;
bool mip_;
int current_solution_index_;
MPCallback* callback_ = nullptr;
@@ -541,7 +545,7 @@ int GUROBI_STDCALL CallbackImpl(GRBmodel* model,
} // namespace
void GurobiInterface::CheckedGurobiCall(int err) const {
::operations_research::CheckedGurobiCall(err, env_);
::operations_research::CheckedGurobiCall(err, global_env_);
}
void GurobiInterface::SetIntAttr(const char* name, int value) {
@@ -608,11 +612,11 @@ char GurobiInterface::GetCharAttrElement(const char* name, int index) const {
GurobiInterface::GurobiInterface(MPSolver* const solver, bool mip)
: MPSolverInterface(solver),
model_(nullptr),
env_(nullptr),
global_env_(nullptr),
mip_(mip),
current_solution_index_(0) {
env_ = GetGurobiEnv().value();
CheckedGurobiCall(GRBnewmodel(env_, &model_, solver_->name_.c_str(),
global_env_ = GetGurobiEnv().value();
CheckedGurobiCall(GRBnewmodel(global_env_, &model_, solver_->name_.c_str(),
0, // numvars
nullptr, // obj
nullptr, // lb
@@ -620,13 +624,13 @@ GurobiInterface::GurobiInterface(MPSolver* const solver, bool mip)
nullptr, // vtype
nullptr)); // varnanes
SetIntAttr(GRB_INT_ATTR_MODELSENSE, maximize_ ? GRB_MAXIMIZE : GRB_MINIMIZE);
CheckedGurobiCall(GRBsetintparam(env_, GRB_INT_PAR_THREADS,
CheckedGurobiCall(GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_THREADS,
absl::GetFlag(FLAGS_num_gurobi_threads)));
}
GurobiInterface::~GurobiInterface() {
CheckedGurobiCall(GRBfreemodel(model_));
GRBfreeenv(env_);
GRBfreeenv(global_env_);
}
// ------ Model modifications and extraction -----
@@ -636,7 +640,7 @@ void GurobiInterface::Reset() {
const absl::MutexLock lock(&hold_interruptions_mutex_);
GRBmodel* old_model = model_;
CheckedGurobiCall(GRBnewmodel(env_, &model_, solver_->name_.c_str(),
CheckedGurobiCall(GRBnewmodel(global_env_, &model_, solver_->name_.c_str(),
0, // numvars
nullptr, // obj
nullptr, // lb
@@ -1210,7 +1214,7 @@ MPSolver::ResultStatus GurobiInterface::Solve(const MPSolverParameters& param) {
CheckedGurobiCall(GRBsetcallbackfunc(model_, nullptr, nullptr));
} else {
gurobi_context = std::make_unique<GurobiMPCallbackContext>(
env_, &mp_var_to_gurobi_var_, num_gurobi_vars_,
global_env_, &mp_var_to_gurobi_var_, num_gurobi_vars_,
callback_->might_add_cuts(), callback_->might_add_lazy_constraints());
mp_callback_with_context.context = gurobi_context.get();
mp_callback_with_context.callback = callback_;
@@ -1235,7 +1239,7 @@ MPSolver::ResultStatus GurobiInterface::Solve(const MPSolverParameters& param) {
const int status = GRBoptimize(model_);
if (status) {
VLOG(1) << "Failed to optimize MIP." << GRBgeterrormsg(env_);
VLOG(1) << "Failed to optimize MIP." << GRBgeterrormsg(global_env_);
} else {
VLOG(1) << absl::StrFormat("Solved in %s.",
absl::FormatDuration(timer.GetDuration()));
@@ -1277,7 +1281,7 @@ MPSolver::ResultStatus GurobiInterface::Solve(const MPSolverParameters& param) {
GRBgetdblattr(model_, GRB_DBL_ATTR_OBJBOUND, &best_objective_bound_);
LOG_IF(WARNING, error != 0)
<< "Best objective bound is not available, error=" << error
<< ", message=" << GRBgeterrormsg(env_);
<< ", message=" << GRBgeterrormsg(global_env_);
VLOG(1) << "best bound = " << best_objective_bound_;
}
@@ -1336,7 +1340,7 @@ std::optional<MPSolutionResponse> GurobiInterface::DirectlySolveProto(
// Here we reuse the Gurobi environment to support single-use license that
// forbids creating a second environment if one already exists.
const auto status_or = GurobiSolveProto(request, env_);
const auto status_or = GurobiSolveProto(request, global_env_);
if (status_or.ok()) return status_or.value();
// Special case: if something is not implemented yet, fall back to solving
// through MPSolver.
@@ -1392,7 +1396,7 @@ void GurobiInterface::Write(const std::string& filename) {
VLOG(1) << "Writing Gurobi model file \"" << filename << "\".";
const int status = GRBwrite(model_, filename.c_str());
if (status) {
LOG(WARNING) << "Failed to write MIP." << GRBgeterrormsg(env_);
LOG(WARNING) << "Failed to write MIP." << GRBgeterrormsg(global_env_);
}
}

35
ortools/math_opt/BUILD.bazel
View File

@@ -11,20 +11,11 @@
# See the License for the specific language governing permissions and
# limitations under the License.
load("@rules_python//python:proto.bzl", "py_proto_library")
load("@rules_cc//cc:defs.bzl", "cc_proto_library")
load("@rules_python//python:proto.bzl", "py_proto_library")
package(default_visibility = ["//visibility:public"])
# Features that are new or under construction have restricted access, contact
# math-opt-dev@ if you want try using these with submitted code.
package_group(
name = "math_opt_allow_list",
packages = [
"//ortools/...",
],
)
proto_library(
name = "callback_proto",
srcs = ["callback.proto"],
@@ -215,3 +206,27 @@ py_proto_library(
name = "infeasible_subsystem_py_pb2",
deps = [":infeasible_subsystem_proto"],
)
proto_library(
name = "rpc_proto",
srcs = ["rpc.proto"],
deps = [
":callback_proto",
":infeasible_subsystem_proto",
":model_parameters_proto",
":model_proto",
":model_update_proto",
":parameters_proto",
":result_proto",
],
)
py_proto_library(
name = "rpc_py_pb2",
deps = [":rpc_proto"],
)
cc_proto_library(
name = "rpc_cc_proto",
deps = [":rpc_proto"],
)

4
ortools/math_opt/core/python/BUILD.bazel
View File

@@ -20,7 +20,7 @@ pybind_extension(
srcs = ["solver.cc"],
deps =
select({
"//ortools/linear_solver:use_scip": ["//ortools/math_opt/solvers:gscip_solver"],
"//ortools/linear_solver:use_cp_sat": ["//ortools/math_opt/solvers:cp_sat_solver"],
"//conditions:default": [],
}) +
select({
@@ -32,7 +32,7 @@ pybind_extension(
"//conditions:default": [],
}) +
select({
"//ortools/linear_solver:use_cp_sat": ["//ortools/math_opt/solvers:cp_sat_solver"],
"//ortools/linear_solver:use_scip": ["//ortools/math_opt/solvers:gscip_solver"],
"//conditions:default": [],
}) + [
"//ortools/math_opt:result_cc_proto",

8
ortools/math_opt/cpp/solve_result.cc
View File

@@ -325,13 +325,13 @@ bool Termination::limit_reached() const {
reason == TerminationReason::kNoSolutionFound;
}
absl::Status Termination::ReasonIs(const TerminationReason reason) const {
absl::Status Termination::EnsureReasonIs(TerminationReason reason) const {
if (this->reason == reason) return absl::OkStatus();
return util::InternalErrorBuilder()
<< "expected termination reason '" << reason << "' but got " << *this;
}
absl::Status Termination::ReasonIsAnyOf(
absl::Status Termination::EnsureReasonIsAnyOf(
std::initializer_list<TerminationReason> reasons) const {
for (const TerminationReason reason : reasons) {
if (this->reason == reason) return absl::OkStatus();
@@ -348,11 +348,11 @@ absl::Status Termination::ReasonIsAnyOf(
}
absl::Status Termination::IsOptimal() const {
return ReasonIs(TerminationReason::kOptimal);
return EnsureReasonIs(TerminationReason::kOptimal);
}
absl::Status Termination::IsOptimalOrFeasible() const {
return ReasonIsAnyOf(
return EnsureReasonIsAnyOf(
{TerminationReason::kOptimal, TerminationReason::kFeasible});
}

4
ortools/math_opt/cpp/solve_result.h
View File

@@ -334,11 +334,11 @@ struct Termination {
// Returns an OkStatus if the reason of this `Termination` is `reason`, or an
// `InternalError` otherwise.
absl::Status ReasonIs(TerminationReason reason) const;
absl::Status EnsureReasonIs(TerminationReason reason) const;
// Returns an OkStatus if the reason of this `Termination` is in `reasons`, or
// an `InternalError` otherwise.
absl::Status ReasonIsAnyOf(
absl::Status EnsureReasonIsAnyOf(
std::initializer_list<TerminationReason> reasons) const;
// Returns termination with reason kOptimal, the provided objective for both

2
ortools/math_opt/python/normalize.py
View File

@@ -20,7 +20,7 @@
"""Utility functions for normalizing proto3 message objects in Python."""
from google3.google.protobuf import duration_pb2
from google.protobuf import duration_pb2
from google.protobuf import descriptor
from google.protobuf import message

61
ortools/math_opt/rpc.proto Normal file
View File

@@ -0,0 +1,61 @@
// Copyright 2010-2022 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.
syntax = "proto3";
package operations_research.math_opt;
import "ortools/math_opt/callback.proto";
import "ortools/math_opt/infeasible_subsystem.proto";
import "ortools/math_opt/model.proto";
import "ortools/math_opt/model_parameters.proto";
import "ortools/math_opt/model_update.proto";
import "ortools/math_opt/parameters.proto";
import "ortools/math_opt/result.proto";
option java_package = "com.google.ortools.mathopt";
option java_multiple_files = true;
// Request for a unary remote solve in MathOpt.
message SolveRequest {
// Solver type to numerically solve the problem. Note that if a solver does
// not support a specific feautre in the model, the optimization procedure
// won't be successful.
SolverTypeProto solver_type = 1;
// A mathematical representation of the optimization problem to solve.
ModelProto model = 2;
SolverInitializerProto initializer = 3;
// Parameters to control a single solve. The enable_output parameter is
// handled specifically. For solvers that support messages callbacks, setting
// it to true will have the server register a message callback. The resulting
// messages will be returned in SolveResponse.messages. For other
// solvers, setting enable_output to true will result in an error.
SolveParametersProto parameters = 4;
// Parameters to control a single solve that are specific to the input model
// (see SolveParametersProto for model independent parameters).
ModelSolveParametersProto model_parameters = 5;
}
// Response for a unary remote solve in MathOpt.
message SolveResponse {
// Description of the output of solving the model in the request.
SolveResultProto result = 1;
// If SolveParametersProto.enable_output has been used, this will contain log
// messages for solvers that support message callbacks.
repeated string messages = 2;
}

3
ortools/math_opt/solvers/BUILD.bazel
View File

@@ -11,8 +11,8 @@
# See the License for the specific language governing permissions and
# limitations under the License.
load("@rules_python//python:proto.bzl", "py_proto_library")
load("@rules_cc//cc:defs.bzl", "cc_proto_library")
load("@rules_python//python:proto.bzl", "py_proto_library")
package(default_visibility = ["//ortools/math_opt:__subpackages__"])
@@ -27,6 +27,7 @@ cc_library(
":gscip_solver_callback",
":message_callback_data",
"//ortools/base:cleanup",
"//ortools/base:linked_hash_map",
"//ortools/base:map_util",
"//ortools/base:protoutil",
"//ortools/base:status_macros",

11
ortools/math_opt/solvers/glop_solver.cc
View File

@@ -353,15 +353,10 @@ absl::StatusOr<glop::GlopParameters> GlopSolver::MergeSolveParameters(
case LP_ALGORITHM_DUAL_SIMPLEX:
result.set_use_dual_simplex(true);
break;
case LP_ALGORITHM_BARRIER:
warnings.emplace_back(
"GLOP does not support 'LP_ALGORITHM_BARRIER' value for "
"'lp_algorithm' parameter.");
break;
default:
LOG(FATAL) << "LPAlgorithm: "
<< ProtoEnumToString(solve_parameters.lp_algorithm())
<< " unknown, error setting GLOP parameters";
warnings.emplace_back(absl::StrCat(
"GLOP does not support the 'lp_algorithm' parameter value: ",
ProtoEnumToString(solve_parameters.lp_algorithm())));
}
}
if (!result.has_use_scaling() && !result.has_scaling_method() &&

72
ortools/math_opt/solvers/gscip_solver.cc
View File

@@ -39,6 +39,7 @@
#include "absl/types/span.h"
#include "google/protobuf/repeated_ptr_field.h"
#include "ortools/base/cleanup.h"
#include "ortools/base/linked_hash_map.h"
#include "ortools/base/logging.h"
#include "ortools/base/map_util.h"
#include "ortools/base/protoutil.h"
@@ -246,60 +247,17 @@ inline GScipVarType GScipVarTypeFromIsInteger(const bool is_integer) {
return is_integer ? GScipVarType::kInteger : GScipVarType::kContinuous;
}
// Used to delay the evaluation of a costly computation until the first time it
// is actually needed.
//
// The typical use is when we have two independent branches that need the same
// data but we don't want to compute these data if we don't enter any of those
// branches.
//
// Usage:
// LazyInitialized<Xxx> xxx([&]() {
// return Xxx(...);
// });
//
// if (predicate_1) {
// ...
// f(xxx.GetOrCreate());
// ...
// }
// if (predicate_2) {
// ...
// f(xxx.GetOrCreate());
// ...
// }
template <typename T>
class LazyInitialized {
public:
// Checks that the input initializer is not nullptr.
explicit LazyInitialized(std::function<T()> initializer)
: initializer_(ABSL_DIE_IF_NULL(initializer)) {}
// Returns the value returned by initializer(), calling it the first time.
const T& GetOrCreate() {
if (!value_) {
value_ = initializer_();
}
return *value_;
}
private:
const std::function<T()> initializer_;
std::optional<T> value_;
};
template <typename T>
SparseDoubleVectorProto FillSparseDoubleVector(
const std::vector<int64_t>& ids_in_order,
const absl::flat_hash_map<int64_t, T>& id_map,
const gtl::linked_hash_map<int64_t, T>& id_map,
const absl::flat_hash_map<T, double>& value_map,
const SparseVectorFilterProto& filter) {
SparseVectorFilterPredicate predicate(filter);
SparseDoubleVectorProto result;
for (const int64_t variable_id : ids_in_order) {
const double value = value_map.at(id_map.at(variable_id));
if (predicate.AcceptsAndUpdate(variable_id, value)) {
result.add_ids(variable_id);
for (const auto [mathopt_id, scip_id] : id_map) {
const double value = value_map.at(scip_id);
if (predicate.AcceptsAndUpdate(mathopt_id, value)) {
result.add_ids(mathopt_id);
result.add_values(value);
}
}
@@ -989,15 +947,6 @@ absl::StatusOr<SolveResultProto> GScipSolver::CreateSolveResultProto(
}
};
LazyInitialized<std::vector<int64_t>> sorted_variables([&]() {
std::vector<int64_t> sorted;
sorted.reserve(variables_.size());
for (const auto& entry : variables_) {
sorted.emplace_back(entry.first);
}
std::sort(sorted.begin(), sorted.end());
return sorted;
});
CHECK_EQ(gscip_result.solutions.size(), gscip_result.objective_values.size());
for (int i = 0; i < gscip_result.solutions.size(); ++i) {
// GScip ensures the solutions are returned best objective first.
@@ -1009,14 +958,13 @@ absl::StatusOr<SolveResultProto> GScipSolver::CreateSolveResultProto(
solution->mutable_primal_solution();
primal_solution->set_objective_value(gscip_result.objective_values[i]);
primal_solution->set_feasibility_status(SOLUTION_STATUS_FEASIBLE);
*primal_solution->mutable_variable_values() = FillSparseDoubleVector(
sorted_variables.GetOrCreate(), variables_, gscip_result.solutions[i],
model_parameters.variable_values_filter());
*primal_solution->mutable_variable_values() =
FillSparseDoubleVector(variables_, gscip_result.solutions[i],
model_parameters.variable_values_filter());
}
if (!gscip_result.primal_ray.empty()) {
*solve_result.add_primal_rays()->mutable_variable_values() =
FillSparseDoubleVector(sorted_variables.GetOrCreate(), variables_,
gscip_result.primal_ray,
FillSparseDoubleVector(variables_, gscip_result.primal_ray,
model_parameters.variable_values_filter());
}
ASSIGN_OR_RETURN(*solve_result.mutable_termination(),

4
ortools/math_opt/solvers/gscip_solver.h
View File

@@ -25,6 +25,7 @@
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/types/span.h"
#include "ortools/base/linked_hash_map.h"
#include "ortools/gscip/gscip.h"
#include "ortools/gscip/gscip.pb.h"
#include "ortools/gscip/gscip_event_handler.h"
@@ -195,8 +196,9 @@ class GScipSolver : public SolverInterface {
absl::Status RegisterHandlers();
const std::unique_ptr<GScip> gscip_;
InterruptEventHandler interrupt_event_handler_;
absl::flat_hash_map<int64_t, SCIP_VAR*> variables_;
gtl::linked_hash_map<int64_t, SCIP_VAR*> variables_;
bool has_quadratic_objective_ = false;
absl::flat_hash_map<int64_t, SCIP_CONS*> linear_constraints_;
absl::flat_hash_map<int64_t, SCIP_CONS*> quadratic_constraints_;

5
ortools/math_opt/solvers/gurobi_solver.cc
View File

@@ -891,6 +891,7 @@ absl::StatusOr<SolveResultProto> GurobiSolver::ExtractSolveResultProto(
GetBestDualBound(solution_and_claims.solutions));
solution_claims = solution_and_claims.solution_claims;
// TODO(b/195295177): Add tests for rays in unbounded MIPs
RETURN_IF_ERROR(FillRays(model_parameters, solution_claims, result));
for (SolutionProto& solution : solution_and_claims.solutions) {
@@ -2929,7 +2930,9 @@ absl::StatusOr<SolveResultProto> GurobiSolver::Solve(
RETURN_IF_ERROR(
UpdateInt32ListAttribute(model_parameters.branching_priorities(),
GRB_INT_ATTR_BRANCHPRIORITY, variables_map_));
RETURN_IF_ERROR(SetMultiObjectiveTolerances(model_parameters));
if (is_multi_objective_mode()) {
RETURN_IF_ERROR(SetMultiObjectiveTolerances(model_parameters));
}
// Here we register the callback when we either have a user callback or a
// local interrupter. The rationale for doing so when we have only an

2
ortools/math_opt/solvers/gurobi_solver.h
View File

@@ -379,6 +379,8 @@ class GurobiSolver : public SolverInterface {
// `.constraint_index` is set and refers to the Gurobi linear constraint index
// for a slack constraint just added to the model such that:
// `expression` == `.variable_index`.
// TODO(b/267310257): Use this for linear constraint slacks, and maybe move it
// up the stack to a bridge.
absl::StatusOr<VariableEqualToExpression>
CreateSlackVariableEqualToExpression(const LinearExpressionProto& expression);

24
ortools/packing/README.md
View File

@@ -13,3 +13,27 @@ constraint and has the maximum value. There can be only one container. More
generally, the knapsack problem models resource allocation.
More information on [Wikipedia](https://en.wikipedia.org/wiki/Knapsack_problem).
Multi-dimensional knapsack solver:
* [`algorithms/knapsack_solver.h`](../algorithms/knapsack_solver.h)
## Bin packing
In an instance of the bin-packing problem, the goal is to use the minimum amount
of bins to fit a series of items of varying size. Vector bin packing (VBP)
considers several dimensions for item size and bin capacity.
More information on
[Wikipedia](https://en.wikipedia.org/wiki/Bin_packing_problem).
Vector bin packing:
* Problem description: [`vector_bin_packing.proto`](vector_bin_packing.proto)
* Solver: [`arc_flow_solver.h`](arc_flow_solver.h)
* VBP parser: [`vector_bin_packing_parser.h`](vector_bin_packing_parser.h)
2D bin packing:
* Problem description:
[`multiple_dimensions_bin_packing.proto`](multiple_dimensions_bin_packing.proto)