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

[CP-SAT] more presolve

Browse Source
This commit is contained in:
Laurent Perron
2025-04-25 19:02:39 +02:00
parent a866dd1b6a
commit 9845bef9fa
11 changed files with 707 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
ortools/sat/BUILD.bazel
View File

@@ -176,11 +176,12 @@ cc_library(
deps = [
":cp_model_cc_proto",
":cp_model_utils",
"//ortools/base:stl_util",
"//ortools/util:bitset",
"//ortools/util:sorted_interval_list",
"@abseil-cpp//absl/algorithm:container",
"@abseil-cpp//absl/container:btree",
"@abseil-cpp//absl/log:check",
"@abseil-cpp//absl/types:span",
],
)
@@ -1718,6 +1719,7 @@ cc_test(
cc_library(
name = "integer_base",
srcs = ["integer_base.cc"],
hdrs = ["integer_base.h"],
deps = [
":sat_base",
@@ -1726,12 +1728,23 @@ cc_library(
"//ortools/util:saturated_arithmetic",
"//ortools/util:sorted_interval_list",
"//ortools/util:strong_integers",
"@abseil-cpp//absl/container:flat_hash_map",
"@abseil-cpp//absl/log:check",
"@abseil-cpp//absl/strings",
"@abseil-cpp//absl/types:span",
],
)
cc_test(
name = "integer_base_test",
size = "small",
srcs = ["integer_base_test.cc"],
deps = [
":integer_base",
"//ortools/base:gmock_main",
],
)
cc_library(
name = "integer",
srcs = ["integer.cc"],

7
ortools/sat/cp_model_checker.cc
View File

@@ -2003,6 +2003,12 @@ bool SolutionIsFeasible(const CpModelProto& model,
VLOG(2) << "Checker inner objective = " << inner_objective;
VLOG(2) << "Checker scaled objective = " << scaled_objective;
}
return true;
}
bool SolutionCanBeOptimal(const CpModelProto& model,
absl::Span<const int64_t> variable_values) {
if (absl::GetFlag(FLAGS_cp_model_check_dependent_variables)) {
const VariableRelationships relationships =
ComputeVariableRelationships(model);
@@ -2015,7 +2021,6 @@ bool SolutionIsFeasible(const CpModelProto& model,
&all_variables));
CHECK(absl::MakeSpan(all_variables) == variable_values);
}
return true;
}

4
ortools/sat/cp_model_checker.h
View File

@@ -68,6 +68,10 @@ bool SolutionIsFeasible(const CpModelProto& model,
const CpModelProto* mapping_proto = nullptr,
const std::vector<int>* postsolve_mapping = nullptr);
// Verifies some invariants that any optimal solution must satisfy.
bool SolutionCanBeOptimal(const CpModelProto& model,
absl::Span<const int64_t> variable_values);
// Checks a single constraint for feasibility.
// This has some overhead, and should only be used for debugging.
// The full model is needed for scheduling constraints that refers to intervals.

132
ortools/sat/cp_model_presolve.cc
View File

@@ -98,6 +98,18 @@ namespace sat {
namespace {
LinearExpression2 GetLinearExpression2FromProto(int a, int64_t coeff_a, int b,
int64_t coeff_b) {
LinearExpression2 result;
DCHECK(RefIsPositive(a));
DCHECK(RefIsPositive(b));
result.vars[0] = IntegerVariable(2 * a);
result.vars[1] = IntegerVariable(2 * b);
result.coeffs[0] = IntegerValue(coeff_a);
result.coeffs[1] = IntegerValue(coeff_b);
return result;
}
// TODO(user): Just make sure this invariant is enforced in all our linear
// constraint after copy, and simplify the code!
bool LinearConstraintIsClean(const LinearConstraintProto& linear) {
@@ -1102,8 +1114,6 @@ bool CpModelPresolver::PresolveLinMax(int c, ConstraintProto* ct) {
unique_var_is_small_enough = context_->DomainOf(unique_var).Size() <= 1000;
}
// This is a test.12y
bool changed;
if (is_one_var_affine_max && unique_var_is_small_enough) {
changed = PropagateAndReduceAffineMax(ct);
@@ -1124,6 +1134,54 @@ bool CpModelPresolver::PresolveLinMax(int c, ConstraintProto* ct) {
return MarkConstraintAsFalse(ct);
}
// Try to reduce lin_max using known relation.
if (ct->lin_max().exprs().size() < 10) {
const int num_exprs = ct->lin_max().exprs().size();
bool simplified = false;
std::vector<bool> can_be_removed(num_exprs, false);
for (int i = 0; i < num_exprs; ++i) {
if (ct->lin_max().exprs(i).vars().size() != 1) continue;
for (int j = 0; j < num_exprs; ++j) {
if (i == j) continue;
if (can_be_removed[j]) continue;
// Note that we skip constant expressions as this should already be
// handled when we compute the domain of each expression and remove
// the ones that are smaller than the target.
if (ct->lin_max().exprs(j).vars().size() != 1) continue;
// Do we know if expr(i) <= expr(j) ?
const LinearExpression2 expr2 = GetLinearExpression2FromProto(
ct->lin_max().exprs(i).vars(0), ct->lin_max().exprs(i).coeffs(0),
ct->lin_max().exprs(j).vars(0), -ct->lin_max().exprs(j).coeffs(0));
const IntegerValue lb = kMinIntegerValue;
const IntegerValue ub(ct->lin_max().exprs(j).offset() -
ct->lin_max().exprs(i).offset());
const RelationStatus status = known_linear2_.GetStatus(expr2, lb, ub);
if (status == RelationStatus::IS_TRUE) {
simplified = true;
can_be_removed[i] = true;
break;
}
}
}
if (simplified) {
context_->UpdateRuleStats(
"lin_max: removed expression smaller than others");
int new_size = 0;
for (int i = 0; i < num_exprs; ++i) {
if (can_be_removed[i]) continue;
*ct->mutable_lin_max()->mutable_exprs(new_size++) =
ct->lin_max().exprs(i);
}
google::protobuf::util::Truncate(ct->mutable_lin_max()->mutable_exprs(),
new_size);
context_->UpdateConstraintVariableUsage(c);
}
}
// If only one is left, we can convert to an equality. Note that we create a
// new constraint otherwise it might not be processed again.
if (ct->lin_max().exprs().size() == 1) {
@@ -2783,6 +2841,36 @@ bool CpModelPresolver::PresolveLinearOfSizeTwo(ConstraintProto* ct) {
const int64_t coeff1 = arg.coeffs(0);
const int64_t coeff2 = arg.coeffs(1);
// Starts by updating our hash map of known relation.
{
const LinearExpression2 expr2 =
GetLinearExpression2FromProto(var1, coeff1, var2, coeff2);
const IntegerValue lb(arg.domain(0));
const IntegerValue ub(arg.domain(arg.domain().size() - 1));
const RelationStatus status = known_linear2_.GetStatus(expr2, lb, ub);
if (status == RelationStatus::IS_TRUE) {
// Note that we don't track what constraint implied the relation, so we
// cannot remove this constraint even if the relation is already known.
//
// However since we only add it if the relation is not
// enforced, this should be correct.
//
// Tricky: If the constraint domain is not simple, we cannot really deduce
// anything.
if (!ct->enforcement_literal().empty() &&
ct->linear().domain().size() == 2) {
context_->UpdateRuleStats("linear2: already known enforced relation");
return RemoveConstraint(ct);
}
} else if (status == RelationStatus::IS_FALSE) {
context_->UpdateRuleStats("linear2: infeasible relation");
return MarkConstraintAsFalse(ct);
} else if (ct->enforcement_literal().empty()) {
known_linear2_.Add(expr2, lb, ub);
}
}
// If it is not an equality, we only presolve the constraint if one of
// the variable is Boolean. Note that if both are Boolean, then a similar
// reduction is done by PresolveLinearOnBooleans(). If we have an equality,
@@ -9931,6 +10019,9 @@ void CpModelPresolver::DetectDuplicateConstraints() {
context_->WriteObjectiveToProto();
}
// If we detect duplicate intervals, we will remap constraints using them.
std::vector<int> interval_mapping;
// Remove duplicate constraints.
// Note that at this point the objective in the proto should be up to date.
//
@@ -9948,6 +10039,13 @@ void CpModelPresolver::DetectDuplicateConstraints() {
? kObjectiveConstraint
: context_->working_model->constraints(rep).constraint_case();
if (type == ConstraintProto::kInterval) {
interval_mapping.resize(context_->working_model->constraints().size(),
-1);
CHECK_EQ(interval_mapping[rep], -1);
interval_mapping[dup] = rep;
}
// For linear constraint, we merge their rhs since it was ignored in the
// FindDuplicateConstraints() call.
if (type == ConstraintProto::kLinear) {
@@ -9998,10 +10096,30 @@ void CpModelPresolver::DetectDuplicateConstraints() {
context_->ReadObjectiveFromProto();
}
}
// Remove the duplicate constraint.
context_->working_model->mutable_constraints(dup)->Clear();
context_->UpdateConstraintVariableUsage(dup);
context_->UpdateRuleStats("duplicate: removed constraint");
}
if (!interval_mapping.empty()) {
context_->UpdateRuleStats("duplicate: remapped duplicate intervals");
const int num_constraints = context_->working_model->constraints().size();
for (int c = 0; c < num_constraints; ++c) {
bool changed = false;
ApplyToAllIntervalIndices(
[&interval_mapping, &changed](int* ref) {
const int new_ref = interval_mapping[*ref];
if (new_ref != -1) {
changed = true;
*ref = new_ref;
}
},
context_->working_model->mutable_constraints(c));
if (changed) context_->UpdateConstraintVariableUsage(c);
}
}
}
void CpModelPresolver::DetectDuplicateConstraintsWithDifferentEnforcements(
@@ -10031,6 +10149,12 @@ void CpModelPresolver::DetectDuplicateConstraintsWithDifferentEnforcements(
auto* dup_ct = context_->working_model->mutable_constraints(dup);
auto* rep_ct = context_->working_model->mutable_constraints(rep);
if (dup_ct->constraint_case() == ConstraintProto::kInterval) {
context_->UpdateRuleStats(
"TODO interval: same interval with different enforcement?");
continue;
}
// Make sure our enforcement list are up to date: nothing fixed and that
// its uses the literal representatives.
if (PresolveEnforcementLiteral(dup_ct)) {
@@ -14202,10 +14326,6 @@ std::vector<std::pair<int, int>> FindDuplicateConstraints(
const auto type = model_proto.constraints(c).constraint_case();
if (type == ConstraintProto::CONSTRAINT_NOT_SET) continue;
// TODO(user): we could delete duplicate identical interval, but we need
// to make sure reference to them are updated.
if (type == ConstraintProto::kInterval) continue;
// Nothing we will presolve in this case.
if (ignore_enforcement && type == ConstraintProto::kBoolAnd) continue;

7
ortools/sat/cp_model_presolve.h
View File

@@ -388,6 +388,13 @@ class CpModelPresolver {
MaxBoundedSubsetSum ub_feasible_;
MaxBoundedSubsetSum ub_infeasible_;
// We have an hash-map of know relation between two variables.
// In particular, this will include all known precedences a <= b.
//
// We reuse an IntegerVariable/IntegerValue based class via
// GetLinearExpression2FromProto() only visible in the .cc.
BestBinaryRelationBounds known_linear2_;
struct IntervalConstraintEq {
const CpModelProto* working_model;
bool operator()(int a, int b) const;

28
ortools/sat/cp_model_solver.cc
View File

@@ -483,15 +483,15 @@ std::string CpModelStats(const CpModelProto& model_proto) {
}
for (const DecisionStrategyProto& strategy : model_proto.search_strategy()) {
absl::StrAppend(
&result, "Search strategy: on ",
strategy.exprs().size() + strategy.variables().size(), " variables, ",
ProtoEnumToString<DecisionStrategyProto::VariableSelectionStrategy>(
strategy.variable_selection_strategy()),
", ",
ProtoEnumToString<DecisionStrategyProto::DomainReductionStrategy>(
strategy.domain_reduction_strategy()),
"\n");
absl::StrAppend(&result, "Search strategy: on ",
strategy.exprs().size() + strategy.variables().size(),
" variables, ",
DecisionStrategyProto::VariableSelectionStrategy_Name(
strategy.variable_selection_strategy()),
", ",
DecisionStrategyProto::DomainReductionStrategy_Name(
strategy.domain_reduction_strategy()),
"\n");
}
auto count_variables_by_type =
@@ -688,8 +688,8 @@ std::string CpSolverResponseStats(const CpSolverResponse& response,
bool has_objective) {
std::string result;
absl::StrAppend(&result, "CpSolverResponse summary:");
absl::StrAppend(&result, "\nstatus: ",
ProtoEnumToString<CpSolverStatus>(response.status()));
absl::StrAppend(&result,
"\nstatus: ", CpSolverStatus_Name(response.status()));
if (has_objective && response.status() != CpSolverStatus::INFEASIBLE) {
absl::StrAppendFormat(&result, "\nobjective: %.16g",
@@ -1710,7 +1710,7 @@ class LnsSolver : public SubSolver {
" [d:", absl::StrFormat("%0.2e", data.difficulty), ", id:", task_id,
", dtime:", absl::StrFormat("%0.2f", data.deterministic_time), "/",
data.deterministic_limit,
", status:", ProtoEnumToString<CpSolverStatus>(data.status),
", status:", CpSolverStatus_Name(data.status),
", #calls:", generator_->num_calls(),
", p:", fully_solved_proportion, "]");
}
@@ -2796,6 +2796,10 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) {
solution_is_feasible =
SolutionIsFeasible(model_proto, response.solution());
}
if (solution_is_feasible && response.status() == CpSolverStatus::OPTIMAL) {
solution_is_feasible =
SolutionCanBeOptimal(model_proto, response.solution());
}
// We dump the response when infeasible, this might help debugging.
if (!solution_is_feasible) {

126
ortools/sat/integer_base.cc Normal file
View File

@@ -0,0 +1,126 @@
// 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/sat/integer_base.h"
#include <cstdint>
#include <numeric>
#include <utility>
#include "absl/log/check.h"
namespace operations_research::sat {
void LinearExpression2::SimpleCanonicalization() {
if (coeffs[0] == 0) vars[0] = kNoIntegerVariable;
if (coeffs[1] == 0) vars[1] = kNoIntegerVariable;
// Corner case when the underlying variable is the same.
if (vars[0] == vars[1]) {
coeffs[0] += coeffs[1];
coeffs[1] = 0;
vars[1] = kNoIntegerVariable;
if (coeffs[0] == 0) vars[0] = kNoIntegerVariable;
}
// Make sure coeff are positive.
for (int i = 0; i < 2; ++i) {
if (coeffs[i] < 0) {
coeffs[i] = -coeffs[i];
vars[i] = NegationOf(vars[i]);
}
}
// Make sure variable are sorted.
if (vars[0] > vars[1]) {
std::swap(vars[0], vars[1]);
std::swap(coeffs[0], coeffs[1]);
}
}
void LinearExpression2::CanonicalizeAndUpdateBounds(IntegerValue& lb,
IntegerValue& ub) {
// We need to be able to negate without overflow.
CHECK_GE(lb, kMinIntegerValue);
CHECK_LE(ub, kMaxIntegerValue);
SimpleCanonicalization();
if (coeffs[0] == 0 || coeffs[1] == 0) return; // abort.
bool negate = false;
if (coeffs[0] == 0) {
if (coeffs[1] != 0) {
negate = !VariableIsPositive(vars[1]);
}
} else {
negate = !VariableIsPositive(vars[0]);
}
if (negate) {
Negate();
std::swap(lb, ub);
lb = -lb;
ub = -ub;
}
// Do gcd division.
const uint64_t gcd = std::gcd(coeffs[0].value(), coeffs[1].value());
if (gcd > 1) {
coeffs[0] /= gcd;
coeffs[1] /= gcd;
ub = FloorRatio(ub, IntegerValue(gcd));
lb = CeilRatio(lb, IntegerValue(gcd));
}
CHECK(coeffs[0] != 0 || vars[0] == kNoIntegerVariable);
CHECK(coeffs[1] != 0 || vars[1] == kNoIntegerVariable);
}
bool BestBinaryRelationBounds::Add(LinearExpression2 expr, IntegerValue lb,
IntegerValue ub) {
expr.CanonicalizeAndUpdateBounds(lb, ub);
if (expr.coeffs[0] == 0 || expr.coeffs[1] == 0) return false;
auto [it, inserted] = best_bounds_.insert({expr, {lb, ub}});
if (inserted) return true;
const auto [known_lb, known_ub] = it->second;
bool restricted = false;
if (lb > known_lb) {
it->second.first = lb;
restricted = true;
}
if (ub < known_ub) {
it->second.second = ub;
restricted = true;
}
return restricted;
}
RelationStatus BestBinaryRelationBounds::GetStatus(LinearExpression2 expr,
IntegerValue lb,
IntegerValue ub) {
expr.CanonicalizeAndUpdateBounds(lb, ub);
if (expr.coeffs[0] == 0 || expr.coeffs[1] == 0) {
return RelationStatus::IS_UNKNOWN;
}
const auto it = best_bounds_.find(expr);
if (it != best_bounds_.end()) {
const auto [known_lb, known_ub] = it->second;
if (lb <= known_lb && ub >= known_ub) return RelationStatus::IS_TRUE;
if (lb > known_ub || ub < known_lb) return RelationStatus::IS_FALSE;
}
return RelationStatus::IS_UNKNOWN;
}
} // namespace operations_research::sat

63
ortools/sat/integer_base.h
View File

@@ -19,11 +19,13 @@
#include <algorithm>
#include <cstdint>
#include <limits>
#include <numeric>
#include <ostream>
#include <string>
#include <utility>
#include <vector>
#include "absl/container/flat_hash_map.h"
#include "absl/log/check.h"
#include "absl/strings/str_cat.h"
#include "absl/types/span.h"
@@ -37,6 +39,12 @@
namespace operations_research {
namespace sat {
// Callbacks that will be called when the search goes back to level 0.
// Callbacks should return false if the propagation fails.
struct LevelZeroCallbackHelper {
std::vector<std::function<bool()>> callbacks;
};
// Value type of an integer variable. An integer variable is always bounded
// on both sides, and this type is also used to store the bounds [lb, ub] of the
// range of each integer variable.
@@ -335,6 +343,61 @@ H AbslHashValue(H h, const AffineExpression& e) {
return h;
}
// A linear expression with at most two variables (coeffs can be zero).
// And some utility to canonicalize them.
struct LinearExpression2 {
// Take the negation of this expression.
void Negate() {
vars[0] = NegationOf(vars[0]);
vars[1] = NegationOf(vars[1]);
}
// This will not change any bounds on the LinearExpression2.
// That is we will not potentially Negate() the expression like
// CanonicalizeAndUpdateBounds() might do.
void SimpleCanonicalization();
// This fully canonicalize this, and update the given bounds accordingly.
void CanonicalizeAndUpdateBounds(IntegerValue& lb, IntegerValue& ub);
bool operator==(const LinearExpression2& o) const {
return vars[0] == o.vars[0] && vars[1] == o.vars[1] &&
coeffs[0] == o.coeffs[0] && coeffs[1] == o.coeffs[1];
}
IntegerValue coeffs[2];
IntegerVariable vars[2];
};
template <typename H>
H AbslHashValue(H h, const LinearExpression2& e) {
h = H::combine(std::move(h), e.vars[0]);
h = H::combine(std::move(h), e.vars[1]);
h = H::combine(std::move(h), e.coeffs[0]);
h = H::combine(std::move(h), e.coeffs[1]);
return h;
}
// Note that we only care about binary relation, not just simple variable bound.
enum class RelationStatus { IS_TRUE, IS_FALSE, IS_UNKNOWN };
class BestBinaryRelationBounds {
public:
// Register the fact that expr \in [lb, ub] is true.
//
// Returns true if this fact is new, that is if the bounds are tighter than
// the current ones.
bool Add(LinearExpression2 expr, IntegerValue lb, IntegerValue ub);
// Returns the known status of expr <= bound.
RelationStatus GetStatus(LinearExpression2 expr, IntegerValue lb,
IntegerValue ub);
private:
// The best bound on the given "canonicalized" expression.
absl::flat_hash_map<LinearExpression2, std::pair<IntegerValue, IntegerValue>>
best_bounds_;
};
// A model singleton that holds the root level integer variable domains.
// we just store a single domain for both var and its negation.
struct IntegerDomains

67
ortools/sat/integer_base_test.cc Normal file
View File

@@ -0,0 +1,67 @@
// 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/sat/integer_base.h"
#include "gtest/gtest.h"
namespace operations_research::sat {
namespace {
TEST(CanonicalizeAffinePrecedenceTest, Basic) {
LinearExpression2 expr;
expr.vars[0] = IntegerVariable(0);
expr.vars[1] = IntegerVariable(2);
expr.coeffs[0] = IntegerValue(4);
expr.coeffs[1] = IntegerValue(2);
IntegerValue lb(0);
IntegerValue ub(11);
expr.CanonicalizeAndUpdateBounds(lb, ub);
EXPECT_EQ(expr.vars[0], IntegerVariable(0));
EXPECT_EQ(expr.vars[1], IntegerVariable(2));
EXPECT_EQ(expr.coeffs[0], IntegerValue(2));
EXPECT_EQ(expr.coeffs[1], IntegerValue(1));
EXPECT_EQ(lb, 0);
EXPECT_EQ(ub, 5);
}
TEST(BestBinaryRelationBoundsTest, Basic) {
LinearExpression2 expr;
expr.vars[0] = IntegerVariable(0);
expr.vars[1] = IntegerVariable(2);
expr.coeffs[0] = IntegerValue(1);
expr.coeffs[1] = IntegerValue(-1);
BestBinaryRelationBounds best_bounds;
EXPECT_TRUE(best_bounds.Add(expr, IntegerValue(0), IntegerValue(5)));
EXPECT_TRUE(best_bounds.Add(expr, IntegerValue(3), IntegerValue(8)));
EXPECT_TRUE(best_bounds.Add(expr, IntegerValue(-1), IntegerValue(4)));
EXPECT_FALSE(
best_bounds.Add(expr, IntegerValue(3), IntegerValue(4))); // best
EXPECT_EQ(RelationStatus::IS_TRUE,
best_bounds.GetStatus(expr, IntegerValue(-10), IntegerValue(4)));
EXPECT_EQ(RelationStatus::IS_TRUE,
best_bounds.GetStatus(expr, IntegerValue(0), IntegerValue(20)));
EXPECT_EQ(RelationStatus::IS_FALSE,
best_bounds.GetStatus(expr, IntegerValue(5), IntegerValue(20)));
EXPECT_EQ(RelationStatus::IS_FALSE,
best_bounds.GetStatus(expr, IntegerValue(-5), IntegerValue(2)));
EXPECT_EQ(RelationStatus::IS_UNKNOWN,
best_bounds.GetStatus(expr, IntegerValue(-5), IntegerValue(3)));
}
} // namespace
} // namespace operations_research::sat

6
ortools/sat/integer_search.h
View File

@@ -119,12 +119,6 @@ struct SearchHeuristics {
// given model.
void ConfigureSearchHeuristics(Model* model);
// Callbacks that will be called when the search goes back to level 0.
// Callbacks should return false if the propagation fails.
struct LevelZeroCallbackHelper {
std::vector<std::function<bool()>> callbacks;
};
// Resets the solver to the given assumptions before calling
// SolveIntegerProblem().
SatSolver::Status ResetAndSolveIntegerProblem(

339
ortools/sat/primary_variables.cc
View File

@@ -15,6 +15,8 @@
#include <algorithm>
#include <cstdint>
#include <cstdlib>
#include <deque>
#include <limits>
#include <tuple>
#include <vector>
@@ -23,10 +25,11 @@
#include "absl/container/btree_map.h"
#include "absl/container/btree_set.h"
#include "absl/log/check.h"
#include "ortools/base/stl_util.h"
#include "absl/types/span.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/cp_model_utils.h"
#include "ortools/util/bitset.h"
#include "ortools/util/sorted_interval_list.h"
namespace operations_research {
namespace sat {
@@ -47,8 +50,8 @@ namespace sat {
// deduce anything, since for some values of u, z and y, the constraint can be
// simplified to x + 3 = x + 3.
void GetRelationshipForConstraint(const ConstraintProto& ct,
std::vector<int>* deducible_vars,
std::vector<int>* input_vars,
absl::btree_set<int>* deducible_vars,
absl::btree_set<int>* input_vars,
int* preferred_to_deduce) {
deducible_vars->clear();
input_vars->clear();
@@ -63,31 +66,23 @@ void GetRelationshipForConstraint(const ConstraintProto& ct,
}
for (int i = 0; i < ct.linear().vars_size(); ++i) {
if (ct.linear().coeffs(i) == 0) continue;
deducible_vars->push_back(ct.linear().vars(i));
deducible_vars->insert(ct.linear().vars(i));
}
// Let's be defensive if this code get called with a linear that is not
// normalized.
gtl::STLSortAndRemoveDuplicates(deducible_vars);
return;
}
case ConstraintProto::kLinMax: {
// We can deduce only the variables that are only in the target.
absl::btree_set<int> deducible_vars_set;
for (int i = 0; i < ct.lin_max().target().vars_size(); ++i) {
if (ct.lin_max().target().coeffs(i) == 0) continue;
deducible_vars_set.insert(ct.lin_max().target().vars(i));
deducible_vars->insert(ct.lin_max().target().vars(i));
}
for (const auto& expr : ct.lin_max().exprs()) {
for (const int var : expr.vars()) {
input_vars->push_back(var);
input_vars->insert(var);
}
}
gtl::STLSortAndRemoveDuplicates(input_vars);
for (const int var : *input_vars) {
deducible_vars_set.erase(var);
}
for (const int var : deducible_vars_set) {
deducible_vars->push_back(var);
deducible_vars->erase(var);
}
return;
}
@@ -105,7 +100,7 @@ void GetRelationshipForConstraint(const ConstraintProto& ct,
if (ct.int_prod().target().coeffs(i) == 0) continue;
const int var = ct.int_prod().target().vars(i);
if (variable_appearance_count[var] == 1) {
deducible_vars->push_back(var);
deducible_vars->insert(var);
}
}
for (const auto& expr : ct.int_prod().exprs()) {
@@ -116,13 +111,13 @@ void GetRelationshipForConstraint(const ConstraintProto& ct,
// coefficient is 1, but we risk trying to deduce x from 0 = 0 * x.
// TODO(user): do it when the target domain doesn't include 0,
// but use preferred_to_deduce to prefer the target.
input_vars->push_back(var);
input_vars->insert(var);
}
}
}
for (const auto& [var, count] : variable_appearance_count) {
if (count != 1) {
input_vars->push_back(var);
input_vars->insert(var);
}
}
return;
@@ -132,6 +127,81 @@ void GetRelationshipForConstraint(const ConstraintProto& ct,
}
}
void CreateLinMaxFromLinearsAndObjective(
const CpModelProto& model, int var_for_target,
absl::Span<const int> linear_constraint_indexes,
bool var_in_objective_is_negative, ConstraintProto* new_constraint) {
// A variable that is in the objective with a positive coefficient and only
// appears in inequalities will be at the lowest value that is greater or
// equal than the variable domain lower bound and that does not violate any
// bound coming from the inequalities. A similar reasoning applies for a
// variable with a negative coefficient in the objective.
LinearArgumentProto& lin_max = *new_constraint->mutable_lin_max();
lin_max.mutable_target()->add_coeffs(var_in_objective_is_negative ? -1 : 1);
lin_max.mutable_target()->add_vars(var_for_target);
const Domain var_domain =
ReadDomainFromProto(model.variables(var_for_target));
// Add the bound coming from the variable domain.
if (var_in_objective_is_negative) {
lin_max.add_exprs()->set_offset(-var_domain.Max());
} else {
lin_max.add_exprs()->set_offset(var_domain.Min());
}
for (const int c : linear_constraint_indexes) {
const LinearConstraintProto& lin = model.constraints(c).linear();
int64_t coeff = 0;
for (int i = 0; i < lin.coeffs_size(); ++i) {
if (lin.vars(i) == var_for_target) {
coeff = lin.coeffs(i);
break;
}
}
LinearExpressionProto* expr = lin_max.add_exprs();
DCHECK_EQ(lin.domain_size(), 2);
const int64_t expr_min = lin.domain(0);
const int64_t expr_max = lin.domain(1);
if ((coeff < 0) == var_in_objective_is_negative) {
expr->set_offset(expr_min);
} else {
expr->set_offset(-expr_max);
}
const int64_t multiplier =
((coeff < 0) == var_in_objective_is_negative) ? -1 : 1;
for (int i = 0; i < lin.coeffs_size(); ++i) {
if (lin.vars(i) == var_for_target) {
continue;
}
expr->add_vars(lin.vars(i));
expr->add_coeffs(multiplier * lin.coeffs(i));
}
}
}
bool IsObjectiveConstraining(const CpModelProto& model) {
if (!model.has_objective()) {
return false;
}
if (model.objective().domain().empty()) {
return false;
}
if (model.objective().domain_size() > 2) {
return true;
}
int64_t lb = 0;
int64_t ub = 0;
for (int i = 0; i < model.objective().vars_size(); ++i) {
const int var = model.objective().vars(i);
lb += model.objective().coeffs(i) * model.variables(var).domain(0);
ub += model.objective().coeffs(i) *
model.variables(var).domain(model.variables(var).domain_size() - 1);
}
if (model.objective().domain(0) > lb || model.objective().domain(1) < ub) {
return true;
}
return false;
}
VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
VariableRelationships result;
Bitset64<int> var_is_secondary(model.variables_size());
@@ -143,25 +213,113 @@ VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
model.variables_size());
struct ConstraintData {
std::vector<int> deducible_vars;
std::vector<int> input_vars;
// These sets will hold only the "undecided" variables. When a variable
// is marked as primary or secondary, it will be removed.
absl::btree_set<int> deducible_vars;
absl::btree_set<int> input_vars;
int preferred_to_deduce;
// A variable is "undecided" until we decide to mark it as primary or
// secondary.
int num_undecided_vars;
// If a variable participates in the model only via linear inequalities and
// the objective, and *all* the other variables in those inequalities are
// already tagged as primary or secondary, then this variable can be flagged
// as a secondary variable and can be computed as a lin_max of the others.
bool is_linear_inequality;
};
std::vector<ConstraintData> constraint_data(model.constraints_size());
std::vector<std::vector<int>> constraints_per_var(model.variables_size());
Bitset64<int> var_appears_only_in_objective_and_linear(
model.variables_size());
Bitset64<int> var_in_objective_is_negative(model.variables_size());
if (!IsObjectiveConstraining(model)) {
// TODO(user): if we have a constraining objective we can suppose a
// non-constraining one + a linear constraint. But this will allow us to
// find at most one new secondary variable, since all the variables in the
// objective will be connected via this linear constraint.
for (int i = 0; i < model.objective().vars_size(); ++i) {
if (model.objective().coeffs(i) == 0) continue;
var_appears_only_in_objective_and_linear.Set(model.objective().vars(i));
if (model.objective().coeffs(i) < 0) {
var_in_objective_is_negative.Set(model.objective().vars(i));
}
}
}
for (int c = 0; c < model.constraints_size(); ++c) {
ConstraintData& data = constraint_data[c];
const ConstraintProto& ct = model.constraints(c);
data.is_linear_inequality = false;
GetRelationshipForConstraint(ct, &data.deducible_vars, &data.input_vars,
&data.preferred_to_deduce);
data.num_undecided_vars =
data.deducible_vars.size() + data.input_vars.size();
// Now prepare the data for the handling the case of variables that only
// appear in the objective and linear inequalities. There are two cases:
// - either the constraint that is one such linear inequality and we flag it
// as such;
// - if not, we flag all the variables used in this constraint as appearing
// in constraints that are not linear inequalities.
if (ct.constraint_case() == ConstraintProto::kLinear &&
data.deducible_vars.empty() && // Not allowing to fix a secondary var
// directly (i.e., an equality)
ct.enforcement_literal().empty() && ct.linear().domain_size() == 2) {
// This is the kind of inequality we might use for a lin_max
data.is_linear_inequality = true;
for (int i = 0; i < ct.linear().vars_size(); ++i) {
const int var = ct.linear().vars(i);
if (!var_appears_only_in_objective_and_linear.IsSet(var)) {
data.input_vars.insert(var);
continue;
}
if (ct.linear().coeffs(i) == 0) continue;
if (std::abs(ct.linear().coeffs(i)) == 1) {
data.deducible_vars.insert(var);
} else {
data.input_vars.insert(var);
// TODO(user): we can support non-unit coefficients to deduce a
// lin_max from the objective. It will become more difficult, since
// first we will need to compute the lcm of all coefficients (and
// avoid overflow). Then, we will need to build a constraint that will
// be div(target, lin_max(exprs) + lcm - 1, lcm).
var_appears_only_in_objective_and_linear.Set(var, false);
}
}
} else {
// Other kind of constraints, tagged those variables as "used elsewhere".
for (const int var : UsedVariables(ct)) {
var_appears_only_in_objective_and_linear.Set(var, false);
}
}
}
// In the loop above, we lazily set some variables as deducible from linear
// inequalities if they only appeared in the objective and linear inequalities
// when we saw the constraint, but we did not checked how they were used in
// following constraints. Now remove them if it was used in other constraints.
std::vector<int> deducible_vars_to_remove;
for (int c = 0; c < model.constraints_size(); ++c) {
deducible_vars_to_remove.clear();
ConstraintData& data = constraint_data[c];
if (!data.is_linear_inequality) {
continue;
}
for (const int var : data.deducible_vars) {
if (!var_appears_only_in_objective_and_linear.IsSet(var)) {
deducible_vars_to_remove.push_back(var);
data.input_vars.insert(var);
}
}
for (const int var : deducible_vars_to_remove) {
data.deducible_vars.erase(var);
}
if (data.deducible_vars.empty()) {
data.is_linear_inequality = false;
}
}
for (int c = 0; c < constraint_data.size(); ++c) {
ConstraintData& data = constraint_data[c];
if (data.deducible_vars.empty()) {
data.input_vars.clear();
continue;
}
if (data.preferred_to_deduce != -1) {
@@ -178,6 +336,44 @@ VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
}
}
// Variables that cannot be potentially deduced using any constraints are
// primary. Flag them as such and strip them from our data structures.
for (int v = 0; v < model.variables_size(); ++v) {
if (num_times_variable_appears_as_deducible[v] != 0) {
continue;
}
num_times_variable_appears_as_input[v] = 0;
var_is_primary.Set(v);
for (const int c : constraints_per_var[v]) {
ConstraintData& data = constraint_data[c];
data.deducible_vars.erase(v);
data.input_vars.erase(v);
}
constraints_per_var[v].clear();
}
// Now, for variables that only appear in the objective and linear
// inequalities, we count how far we are from being able to deduce their
// value. In practice, we count the number of linear inequalities in which
// this variable appears alongside another variable that we have not decided
// to be primary or secondary yet. When this count reach 0, it means we can
// create a lin_max constraint to deduce its value.
std::vector<int> count_of_unresolved_linear_inequalities_per_var(
model.variables_size());
for (int c = 0; c < model.constraints_size(); ++c) {
ConstraintData& data = constraint_data[c];
if (!data.is_linear_inequality) {
continue;
}
if (data.input_vars.size() + data.deducible_vars.size() > 1) {
for (const int v : data.input_vars) {
count_of_unresolved_linear_inequalities_per_var[v]++;
}
for (const int v : data.deducible_vars) {
count_of_unresolved_linear_inequalities_per_var[v]++;
}
}
}
// Now do a greedy heuristic: we will take each variable in a particular
// order, and if the variable can be deduced from other variables that we have
// already decided to declare as primary or secondary, we will mark it as
@@ -188,11 +384,14 @@ VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
// we reach a fixed point. The heuristic for the order is to try to process
// first the variables that are more "useful" to be marked as primary, so it
// allows us to mark more variables as secondary in the following.
std::vector<int> variables_ordered_by_preference(model.variables_size());
std::deque<int> vars_queue;
for (int v = 0; v < model.variables_size(); ++v) {
variables_ordered_by_preference[v] = v;
if (var_is_primary.IsSet(v) || var_is_secondary.IsSet(v)) {
continue;
}
vars_queue.push_back(v);
}
absl::c_stable_sort(variables_ordered_by_preference, [&](int a, int b) {
absl::c_stable_sort(vars_queue, [&](int a, int b) {
const int a_diff_usage = num_times_variable_appears_as_deducible[a] -
num_times_variable_appears_as_input[a];
const int b_diff_usage = num_times_variable_appears_as_deducible[b] -
@@ -207,14 +406,20 @@ VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
-num_times_variable_appears_as_deducible[b]);
});
std::vector<int> constraints_to_check;
for (int v : variables_ordered_by_preference) {
while (!vars_queue.empty()) {
const int v = vars_queue.front();
vars_queue.pop_front();
if (var_is_secondary.IsSet(v) || var_is_primary.IsSet(v)) {
continue;
}
// first, we will decide whether we should mark `v` as secondary or primary.
// First, we will decide whether we should mark `v` as secondary or primary
// using the equality constraints.
for (const int c : constraints_per_var[v]) {
ConstraintData& data = constraint_data[c];
if (data.num_undecided_vars != 1) {
if (data.is_linear_inequality) {
continue;
}
if (data.deducible_vars.size() + data.input_vars.size() != 1) {
// One of it inputs are neither primary nor secondary yet, we cannot
// deduce the value of `v` using this constraint.
continue;
@@ -222,11 +427,11 @@ VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
// There is a single undecided variable for this constraint. Thus `v` is
// either an input or a deducible variable. Let's check.
if (!absl::c_linear_search(data.input_vars, v)) {
if (data.deducible_vars.empty()) {
// This is a strange case, like `z = lin_max(x, y)`, where `z` and `y`
// are primary (we cannot deduce `x`). Let's just flag this constraint
// as useless from now on.
data.num_undecided_vars = 0;
data.input_vars.clear();
continue;
}
var_is_secondary.Set(v);
@@ -236,30 +441,55 @@ VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
break;
}
// We couldn't deduce the value of `v` from any constraint, so we mark it as
// primary.
// We couldn't deduce the value of `v` from any constraint, check if it only
// appears in linear inequalities.
if (!var_is_secondary.IsSet(v)) {
var_is_primary.Set(v);
if (var_appears_only_in_objective_and_linear.IsSet(v) &&
count_of_unresolved_linear_inequalities_per_var[v] == 0) {
var_is_secondary.Set(v);
result.secondary_variables.push_back(v);
result.dependency_resolution_constraint.emplace_back();
CreateLinMaxFromLinearsAndObjective(
model, v, constraints_per_var[v],
var_in_objective_is_negative.IsSet(v),
&result.dependency_resolution_constraint.back());
// TODO(user): set result.redundant_constraint_indices?
} else {
// We can't deduce the value of `v` from what we have so far, flag it as
// primary.
var_is_primary.Set(v);
}
}
// At this point we have decided to make `v` primary or secondary, so we
// can remove it from the model and maybe lazily deduce some variables.
auto update_constraints_after_var_is_decided = [&](int v) {
for (const int c : constraints_per_var[v]) {
ConstraintData& data = constraint_data[c];
data.num_undecided_vars--;
if (data.num_undecided_vars != 1) {
data.deducible_vars.erase(v);
data.input_vars.erase(v);
if (data.input_vars.size() + data.deducible_vars.size() != 1) {
// Two of the variables for this constraint are neither primary nor
// secondary yet, we cannot deduce the value of anything using this
// constraint.
continue;
}
if (absl::c_all_of(data.deducible_vars, [&var_is_secondary,
&var_is_primary](int v) {
return var_is_secondary.IsSet(v) || var_is_primary.IsSet(v);
})) {
// Same as the test for the case `z = lin_max(x, y)` above.
data.num_undecided_vars = 0;
if (data.is_linear_inequality && data.deducible_vars.size() == 1) {
const int deducible_var = *data.deducible_vars.begin();
count_of_unresolved_linear_inequalities_per_var[deducible_var]--;
if (count_of_unresolved_linear_inequalities_per_var[deducible_var] ==
0) {
// Now we can deduce a new variable from linears, process it ASAP!
vars_queue.push_front(deducible_var);
}
} else {
constraints_to_check.push_back(c);
if (data.deducible_vars.empty()) {
// Same as the test for the case `z = lin_max(x, y)` above.
data.input_vars.clear();
} else {
// This constraint fix a secondary variable, enqueue it!
constraints_to_check.push_back(c);
}
}
}
};
@@ -275,24 +505,11 @@ VariableRelationships ComputeVariableRelationships(const CpModelProto& model) {
const int c = constraints_to_check.back();
constraints_to_check.pop_back();
ConstraintData& data = constraint_data[c];
DCHECK_LE(data.num_undecided_vars, 1);
if (data.num_undecided_vars != 1) {
continue;
}
int single_deducible_var = -1;
for (const int deducible_var : data.deducible_vars) {
if (var_is_secondary.IsSet(deducible_var) ||
var_is_primary.IsSet(deducible_var)) {
continue;
}
DCHECK_EQ(single_deducible_var, -1);
single_deducible_var = deducible_var;
}
if (single_deducible_var == -1) {
// Same as the test for the case `z = lin_max(x, y)` above.
data.num_undecided_vars = 0;
DCHECK_LE(data.deducible_vars.size(), 1);
if (data.deducible_vars.size() != 1) {
continue;
}
const int single_deducible_var = *data.deducible_vars.begin();
var_is_secondary.Set(single_deducible_var);
update_constraints_after_var_is_decided(single_deducible_var);
result.secondary_variables.push_back(single_deducible_var);