[CP-SAT] improve presolve; improve linear subsystem

This commit is contained in:
Laurent Perron
2019-12-16 12:34:56 +01:00
parent bd5c282a51
commit f9f2d7b5b4
22 changed files with 617 additions and 88 deletions

View File

@@ -612,6 +612,28 @@ Constraint CpModelBuilder::AddMinEquality(IntVar target,
return Constraint(proto);
}
void CpModelBuilder::LinearExprToProto(const LinearExpr& expr,
LinearExpressionProto* expr_proto) {
for (const IntVar var : expr.variables()) {
expr_proto->add_vars(GetOrCreateIntegerIndex(var.index_));
}
for (const int64 coeff : expr.coefficients()) {
expr_proto->add_coeffs(coeff);
}
expr_proto->set_offset(expr.constant());
}
Constraint CpModelBuilder::AddLinMinEquality(
const LinearExpr& target, absl::Span<const LinearExpr> exprs) {
ConstraintProto* const proto = cp_model_.add_constraints();
LinearExprToProto(target, proto->mutable_lin_min()->mutable_target());
for (const LinearExpr& expr : exprs) {
LinearExpressionProto* expr_proto = proto->mutable_lin_min()->add_exprs();
LinearExprToProto(expr, expr_proto);
}
return Constraint(proto);
}
Constraint CpModelBuilder::AddMaxEquality(IntVar target,
absl::Span<const IntVar> vars) {
ConstraintProto* const proto = cp_model_.add_constraints();
@@ -622,6 +644,17 @@ Constraint CpModelBuilder::AddMaxEquality(IntVar target,
return Constraint(proto);
}
Constraint CpModelBuilder::AddLinMaxEquality(
const LinearExpr& target, absl::Span<const LinearExpr> exprs) {
ConstraintProto* const proto = cp_model_.add_constraints();
LinearExprToProto(target, proto->mutable_lin_max()->mutable_target());
for (const LinearExpr& expr : exprs) {
LinearExpressionProto* expr_proto = proto->mutable_lin_max()->add_exprs();
LinearExprToProto(expr, expr_proto);
}
return Constraint(proto);
}
Constraint CpModelBuilder::AddDivisionEquality(IntVar target, IntVar numerator,
IntVar denominator) {
ConstraintProto* const proto = cp_model_.add_constraints();

View File

@@ -737,9 +737,17 @@ class CpModelBuilder {
/// Adds target == min(vars).
Constraint AddMinEquality(IntVar target, absl::Span<const IntVar> vars);
/// Adds target == min(exprs).
Constraint AddLinMinEquality(const LinearExpr& target,
absl::Span<const LinearExpr> exprs);
/// Adds target == max(vars).
Constraint AddMaxEquality(IntVar target, absl::Span<const IntVar> vars);
/// Adds target == max(exprs).
Constraint AddLinMaxEquality(const LinearExpr& target,
absl::Span<const LinearExpr> exprs);
/// Adds target = num / denom (integer division rounded towards 0).
Constraint AddDivisionEquality(IntVar target, IntVar numerator,
IntVar denominator);
@@ -812,6 +820,10 @@ class CpModelBuilder {
friend class CumulativeConstraint;
friend class ReservoirConstraint;
// Fills the 'expr_proto' with the linear expression represented by 'expr'.
void LinearExprToProto(const LinearExpr& expr,
LinearExpressionProto* expr_proto);
// Returns a (cached) integer variable index with a constant value.
int IndexFromConstant(int64 value);

View File

@@ -76,6 +76,17 @@ message IntegerArgumentProto {
repeated int32 vars = 2;
}
message LinearExpressionProto {
repeated int32 vars = 1;
repeated int32 coeffs = 2;
int64 offset = 3;
}
message LinearArgumentProto {
LinearExpressionProto target = 1;
repeated LinearExpressionProto exprs = 2;
}
// All variables must take different values.
message AllDifferentConstraintProto {
repeated int32 vars = 1;
@@ -246,7 +257,7 @@ message AutomatonConstraintProto {
repeated int32 vars = 7;
}
// Next id: 27
// Next id: 29
message ConstraintProto {
// For debug/logging only. Can be empty.
string name = 1;
@@ -297,12 +308,22 @@ message ConstraintProto {
// The int_max constraint forces the target to equal the maximum of all
// variables.
// TODO(user): Remove int_max in favor of lin_max.
IntegerArgumentProto int_max = 9;
// The lin_max constraint forces the target to equal the maximum of all
// linear expressions.
LinearArgumentProto lin_max = 27;
// The int_min constraint forces the target to equal the minimum of all
// variables.
// TODO(user): Remove int_min in favor of lin_min.
IntegerArgumentProto int_min = 10;
// The lin_max constraint forces the target to equal the minimum of all
// linear expressions.
LinearArgumentProto lin_min = 28;
// The int_min constraint forces the target to equal the product of all
// variables.
IntegerArgumentProto int_prod = 11;

View File

@@ -194,6 +194,19 @@ std::string ValidateLinearConstraint(const CpModelProto& model,
return "";
}
std::string ValidateLinearExpression(const CpModelProto& model,
const LinearExpressionProto& expr) {
if (expr.coeffs_size() != expr.vars_size()) {
return absl::StrCat("coeffs_size() != vars_size() in linear expression: ",
ProtobufShortDebugString(expr));
}
if (PossibleIntegerOverflow(model, expr)) {
return absl::StrCat("Possible overflow in linear expression: ",
ProtobufShortDebugString(expr));
}
return "";
}
std::string ValidateCircuitConstraint(const CpModelProto& model,
const ConstraintProto& ct) {
const int size = ct.circuit().tails().size();
@@ -387,6 +400,29 @@ std::string ValidateCpModel(const CpModelProto& model) {
}
RETURN_IF_NOT_EMPTY(ValidateLinearConstraint(model, ct));
break;
case ConstraintProto::ConstraintCase::kLinMax: {
const std::string target_error =
ValidateLinearExpression(model, ct.lin_min().target());
if (!target_error.empty()) return target_error;
for (int i = 0; i < ct.lin_max().exprs_size(); ++i) {
const std::string expr_error =
ValidateLinearExpression(model, ct.lin_max().exprs(i));
if (!expr_error.empty()) return expr_error;
}
break;
}
case ConstraintProto::ConstraintCase::kLinMin: {
const std::string target_error =
ValidateLinearExpression(model, ct.lin_min().target());
if (!target_error.empty()) return target_error;
for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
const std::string expr_error =
ValidateLinearExpression(model, ct.lin_min().exprs(i));
if (!expr_error.empty()) return expr_error;
}
break;
}
case ConstraintProto::ConstraintCase::kInterval:
support_enforcement = true;
RETURN_IF_NOT_EMPTY(ValidateIntervalConstraint(model, ct));
@@ -524,6 +560,25 @@ class ConstraintChecker {
return max == actual_max;
}
int64 LinearExpressionValue(const LinearExpressionProto& expr) {
int64 sum = expr.offset();
const int num_variables = expr.vars_size();
for (int i = 0; i < num_variables; ++i) {
sum += Value(expr.vars(i)) * expr.coeffs(i);
}
return sum;
}
bool LinMaxConstraintIsFeasible(const ConstraintProto& ct) {
const int64 max = LinearExpressionValue(ct.lin_max().target());
int64 actual_max = kint64min;
for (int i = 0; i < ct.lin_max().exprs_size(); ++i) {
const int64 expr_value = LinearExpressionValue(ct.lin_max().exprs(i));
actual_max = std::max(actual_max, expr_value);
}
return max == actual_max;
}
bool IntProdConstraintIsFeasible(const ConstraintProto& ct) {
const int64 prod = Value(ct.int_prod().target());
int64 actual_prod = 1;
@@ -552,6 +607,16 @@ class ConstraintChecker {
return min == actual_min;
}
bool LinMinConstraintIsFeasible(const ConstraintProto& ct) {
const int64 min = LinearExpressionValue(ct.lin_min().target());
int64 actual_min = kint64max;
for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
const int64 expr_value = LinearExpressionValue(ct.lin_min().exprs(i));
actual_min = std::min(actual_min, expr_value);
}
return min == actual_min;
}
bool AllDiffConstraintIsFeasible(const ConstraintProto& ct) {
absl::flat_hash_set<int64> values;
for (const int v : ct.all_diff().vars()) {
@@ -949,9 +1014,15 @@ bool SolutionIsFeasible(const CpModelProto& model,
case ConstraintProto::ConstraintCase::kIntMin:
is_feasible = checker.IntMinConstraintIsFeasible(ct);
break;
case ConstraintProto::ConstraintCase::kLinMin:
is_feasible = checker.LinMinConstraintIsFeasible(ct);
break;
case ConstraintProto::ConstraintCase::kIntMax:
is_feasible = checker.IntMaxConstraintIsFeasible(ct);
break;
case ConstraintProto::ConstraintCase::kLinMax:
is_feasible = checker.LinMaxConstraintIsFeasible(ct);
break;
case ConstraintProto::ConstraintCase::kAllDiff:
is_feasible = checker.AllDiffConstraintIsFeasible(ct);
break;

View File

@@ -817,6 +817,33 @@ void ExpandAutomaton(ConstraintProto* ct, PresolveContext* context) {
ct->Clear();
}
// Fills the target as negated ref.
void NegatedLinearExpression(const LinearExpressionProto& ref,
LinearExpressionProto* target) {
for (int i = 0; i < ref.vars_size(); ++i) {
target->add_vars(NegatedRef(ref.vars(i)));
target->add_coeffs(ref.coeffs(i));
}
target->set_offset(-ref.offset());
}
void ExpandLinMax(ConstraintProto* ct, PresolveContext* context) {
ConstraintProto* const lin_min = context->working_model->add_constraints();
for (int i = 0; i < ct->enforcement_literal_size(); ++i) {
lin_min->add_enforcement_literal(ct->enforcement_literal(i));
}
// Target
NegatedLinearExpression(ct->lin_max().target(),
lin_min->mutable_lin_min()->mutable_target());
for (int i = 0; i < ct->lin_max().exprs_size(); ++i) {
LinearExpressionProto* const expr = lin_min->mutable_lin_min()->add_exprs();
NegatedLinearExpression(ct->lin_max().exprs(i), expr);
}
ct->Clear();
}
} // namespace
void ExpandCpModel(PresolveOptions options, PresolveContext* context) {
@@ -838,6 +865,9 @@ void ExpandCpModel(PresolveOptions options, PresolveContext* context) {
case ConstraintProto::ConstraintCase::kIntProd:
ExpandIntProd(ct, context);
break;
case ConstraintProto::ConstraintCase::kLinMax:
ExpandLinMax(ct, context);
break;
case ConstraintProto::ConstraintCase::kElement:
if (options.parameters.expand_element_constraints()) {
ExpandElement(ct, context);

View File

@@ -711,7 +711,9 @@ WeightedRandomRelaxationNeighborhoodGenerator::
case ConstraintProto::kIntDiv:
case ConstraintProto::kIntMod:
case ConstraintProto::kIntMax:
case ConstraintProto::kLinMax:
case ConstraintProto::kIntMin:
case ConstraintProto::kLinMin:
case ConstraintProto::kNoOverlap:
case ConstraintProto::kNoOverlap2D:
constraint_weights_.push_back(2.0);

View File

@@ -19,6 +19,7 @@
#include <set>
#include <string>
#include <utility>
#include <vector>
#include "absl/container/flat_hash_map.h"
#include "absl/container/flat_hash_set.h"
@@ -30,6 +31,7 @@
#include "ortools/sat/all_different.h"
#include "ortools/sat/circuit.h"
#include "ortools/sat/cp_constraints.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/cp_model_utils.h"
#include "ortools/sat/cumulative.h"
#include "ortools/sat/diffn.h"
@@ -1249,6 +1251,29 @@ void LoadIntMinConstraint(const ConstraintProto& ct, Model* m) {
m->Add(IsEqualToMinOf(min, vars));
}
LinearExpression GetExprFromProto(const LinearExpressionProto& expr_proto,
const CpModelMapping& mapping) {
LinearExpression expr;
expr.vars = mapping.Integers(expr_proto.vars());
for (int j = 0; j < expr_proto.coeffs_size(); ++j) {
expr.coeffs.push_back(IntegerValue(expr_proto.coeffs(j)));
}
expr.offset = IntegerValue(expr_proto.offset());
return CanonicalizeExpr(expr);
}
void LoadLinMinConstraint(const ConstraintProto& ct, Model* m) {
auto* mapping = m->GetOrCreate<CpModelMapping>();
const LinearExpression min =
GetExprFromProto(ct.lin_min().target(), *mapping);
std::vector<LinearExpression> exprs;
exprs.reserve(ct.lin_min().exprs_size());
for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
exprs.push_back(GetExprFromProto(ct.lin_min().exprs(i), *mapping));
}
m->Add(IsEqualToMinOf(min, exprs));
}
void LoadIntMaxConstraint(const ConstraintProto& ct, Model* m) {
auto* mapping = m->GetOrCreate<CpModelMapping>();
const IntegerVariable max = mapping->Integer(ct.int_max().target());
@@ -1738,6 +1763,9 @@ bool LoadConstraint(const ConstraintProto& ct, Model* m) {
case ConstraintProto::ConstraintProto::kIntMin:
LoadIntMinConstraint(ct, m);
return true;
case ConstraintProto::ConstraintProto::kLinMin:
LoadLinMinConstraint(ct, m);
return true;
case ConstraintProto::ConstraintProto::kIntMax:
LoadIntMaxConstraint(ct, m);
return true;

View File

@@ -227,6 +227,7 @@ void LoadAllDiffConstraint(const ConstraintProto& ct, Model* m);
void LoadIntProdConstraint(const ConstraintProto& ct, Model* m);
void LoadIntDivConstraint(const ConstraintProto& ct, Model* m);
void LoadIntMinConstraint(const ConstraintProto& ct, Model* m);
void LoadLinMinConstraint(const ConstraintProto& ct, Model* m);
void LoadIntMaxConstraint(const ConstraintProto& ct, Model* m);
void LoadNoOverlapConstraint(const ConstraintProto& ct, Model* m);
void LoadNoOverlap2dConstraint(const ConstraintProto& ct, Model* m);

View File

@@ -104,6 +104,7 @@ class CpModelPresolver {
bool PresolveInterval(int c, ConstraintProto* ct);
bool PresolveIntDiv(ConstraintProto* ct);
bool PresolveIntProd(ConstraintProto* ct);
// TODO(user) : Add presolve rules for LinMin.
bool PresolveIntMin(ConstraintProto* ct);
bool PresolveIntMax(ConstraintProto* ct);
bool PresolveBoolXor(ConstraintProto* ct);

View File

@@ -1968,6 +1968,7 @@ class LnsSolver : public SubSolver {
SatParameters local_params(parameters_);
local_params.set_max_deterministic_time(data.deterministic_limit);
local_params.set_stop_after_first_solution(false);
local_params.set_log_search_progress(false);
if (FLAGS_cp_model_dump_lns) {
const std::string name =

View File

@@ -60,10 +60,24 @@ IndexReferences GetReferencesUsedByConstraint(const ConstraintProto& ct) {
output.variables.push_back(ct.int_max().target());
AddIndices(ct.int_max().vars(), &output.variables);
break;
case ConstraintProto::ConstraintCase::kLinMax: {
AddIndices(ct.lin_max().target().vars(), &output.variables);
for (int i = 0; i < ct.lin_max().exprs_size(); ++i) {
AddIndices(ct.lin_max().exprs(i).vars(), &output.variables);
}
break;
}
case ConstraintProto::ConstraintCase::kIntMin:
output.variables.push_back(ct.int_min().target());
AddIndices(ct.int_min().vars(), &output.variables);
break;
case ConstraintProto::ConstraintCase::kLinMin: {
AddIndices(ct.lin_min().target().vars(), &output.variables);
for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
AddIndices(ct.lin_min().exprs(i).vars(), &output.variables);
}
break;
}
case ConstraintProto::ConstraintCase::kIntProd:
output.variables.push_back(ct.int_prod().target());
AddIndices(ct.int_prod().vars(), &output.variables);
@@ -155,8 +169,12 @@ void ApplyToAllLiteralIndices(const std::function<void(int*)>& f,
break;
case ConstraintProto::ConstraintCase::kIntMax:
break;
case ConstraintProto::ConstraintCase::kLinMax:
break;
case ConstraintProto::ConstraintCase::kIntMin:
break;
case ConstraintProto::ConstraintCase::kLinMin:
break;
case ConstraintProto::ConstraintCase::kIntProd:
break;
case ConstraintProto::ConstraintCase::kLinear:
@@ -218,10 +236,22 @@ void ApplyToAllVariableIndices(const std::function<void(int*)>& f,
APPLY_TO_SINGULAR_FIELD(int_max, target);
APPLY_TO_REPEATED_FIELD(int_max, vars);
break;
case ConstraintProto::ConstraintCase::kLinMax:
APPLY_TO_REPEATED_FIELD(lin_max, target()->mutable_vars);
for (int i = 0; i < ct->lin_max().exprs_size(); ++i) {
APPLY_TO_REPEATED_FIELD(lin_max, exprs(i)->mutable_vars);
}
break;
case ConstraintProto::ConstraintCase::kIntMin:
APPLY_TO_SINGULAR_FIELD(int_min, target);
APPLY_TO_REPEATED_FIELD(int_min, vars);
break;
case ConstraintProto::ConstraintCase::kLinMin:
APPLY_TO_REPEATED_FIELD(lin_min, target()->mutable_vars);
for (int i = 0; i < ct->lin_min().exprs_size(); ++i) {
APPLY_TO_REPEATED_FIELD(lin_min, exprs(i)->mutable_vars);
}
break;
case ConstraintProto::ConstraintCase::kIntProd:
APPLY_TO_SINGULAR_FIELD(int_prod, target);
APPLY_TO_REPEATED_FIELD(int_prod, vars);
@@ -292,8 +322,12 @@ void ApplyToAllIntervalIndices(const std::function<void(int*)>& f,
break;
case ConstraintProto::ConstraintCase::kIntMax:
break;
case ConstraintProto::ConstraintCase::kLinMax:
break;
case ConstraintProto::ConstraintCase::kIntMin:
break;
case ConstraintProto::ConstraintCase::kLinMin:
break;
case ConstraintProto::ConstraintCase::kIntProd:
break;
case ConstraintProto::ConstraintCase::kLinear:
@@ -353,8 +387,12 @@ std::string ConstraintCaseName(
return "kIntMod";
case ConstraintProto::ConstraintCase::kIntMax:
return "kIntMax";
case ConstraintProto::ConstraintCase::kLinMax:
return "kLinMax";
case ConstraintProto::ConstraintCase::kIntMin:
return "kIntMin";
case ConstraintProto::ConstraintCase::kLinMin:
return "kLinMin";
case ConstraintProto::ConstraintCase::kIntProd:
return "kIntProd";
case ConstraintProto::ConstraintCase::kLinear:
@@ -422,8 +460,12 @@ std::vector<int> UsedIntervals(const ConstraintProto& ct) {
break;
case ConstraintProto::ConstraintCase::kIntMax:
break;
case ConstraintProto::ConstraintCase::kLinMax:
break;
case ConstraintProto::ConstraintCase::kIntMin:
break;
case ConstraintProto::ConstraintCase::kLinMin:
break;
case ConstraintProto::ConstraintCase::kIntProd:
break;
case ConstraintProto::ConstraintCase::kLinear:

View File

@@ -791,7 +791,7 @@ void IntegerTrail::AppendRelaxedLinearReason(
for (const IntegerVariable var : vars) {
tmp_indices_.push_back(vars_[var].current_trail_index);
}
RelaxLinearReason(slack, coeffs, &tmp_indices_);
if (slack > 0) RelaxLinearReason(slack, coeffs, &tmp_indices_);
for (const int i : tmp_indices_) {
reason->push_back(IntegerLiteral::GreaterOrEqual(integer_trail_[i].var,
integer_trail_[i].bound));
@@ -810,7 +810,7 @@ void IntegerTrail::RelaxLinearReason(IntegerValue slack,
// We start by filtering *trail_indices:
// - remove all level zero entries.
// - keep the one that cannot be relaxed.
// - move the other one the the relax_heap_ (and creating the heap).
// - move the other one to the relax_heap_ (and creating the heap).
int new_size = 0;
const int size = coeffs.size();
const int num_vars = vars_.size();

View File

@@ -15,9 +15,12 @@
#include <algorithm>
#include <memory>
#include <vector>
#include "absl/container/flat_hash_map.h"
#include "absl/memory/memory.h"
#include "ortools/base/stl_util.h"
#include "ortools/sat/integer.h"
#include "ortools/util/sorted_interval_list.h"
namespace operations_research {
@@ -291,6 +294,146 @@ void MinPropagator::RegisterWith(GenericLiteralWatcher* watcher) {
watcher->WatchUpperBound(min_var_, id);
}
LinearExpression CanonicalizeExpr(const LinearExpression& expr) {
LinearExpression canonical_expr;
canonical_expr.offset = expr.offset;
for (int i = 0; i < expr.vars.size(); ++i) {
if (expr.coeffs[i] < 0) {
canonical_expr.vars.push_back(NegationOf(expr.vars[i]));
canonical_expr.coeffs.push_back(-expr.coeffs[i]);
} else {
canonical_expr.vars.push_back(expr.vars[i]);
canonical_expr.coeffs.push_back(expr.coeffs[i]);
}
}
return canonical_expr;
}
IntegerValue LinExprLowerBound(const LinearExpression& expr,
const IntegerTrail& integer_trail) {
IntegerValue lower_bound = expr.offset;
for (int i = 0; i < expr.vars.size(); ++i) {
DCHECK_GE(expr.coeffs[i], 0) << "The expression is not canonicalized";
lower_bound += expr.coeffs[i] * integer_trail.LowerBound(expr.vars[i]);
}
return lower_bound;
}
IntegerValue LinExprUpperBound(const LinearExpression& expr,
const IntegerTrail& integer_trail) {
IntegerValue upper_bound = expr.offset;
for (int i = 0; i < expr.vars.size(); ++i) {
DCHECK_GE(expr.coeffs[i], 0) << "The expression is not canonicalized";
upper_bound += expr.coeffs[i] * integer_trail.UpperBound(expr.vars[i]);
}
return upper_bound;
}
LinMinPropagator::LinMinPropagator(const std::vector<LinearExpression>& exprs,
IntegerVariable min_var, Model* model)
: exprs_(exprs),
min_var_(min_var),
model_(model),
integer_trail_(model_->GetOrCreate<IntegerTrail>()) {}
bool LinMinPropagator::Propagate() {
ub_constraints_.resize(rev_ub_constraints_size_);
if (exprs_.empty()) return true;
expr_lbs_.clear();
// Count the number of interval that are possible candidate for the min.
// Only the intervals for which lb > current_min_ub cannot.
const IntegerLiteral min_ub_literal =
integer_trail_->UpperBoundAsLiteral(min_var_);
const IntegerValue current_min_ub = integer_trail_->UpperBound(min_var_);
int num_intervals_that_can_be_min = 0;
int last_possible_min_interval = 0;
IntegerValue min_of_linear_expression_lb = kMaxIntegerValue;
for (int i = 0; i < exprs_.size(); ++i) {
const IntegerValue lb = LinExprLowerBound(exprs_[i], *integer_trail_);
expr_lbs_.push_back(lb);
min_of_linear_expression_lb = std::min(min_of_linear_expression_lb, lb);
if (lb <= current_min_ub) {
++num_intervals_that_can_be_min;
last_possible_min_interval = i;
}
}
// Propagation a) lb(min) >= lb(MIN(exprs)) = MIN(lb(expr));
// Conflict will be detected by the fact that the [lb, ub] of the min is
// empty. In case of conflict, we just need the reason for pushing UB + 1.
if (min_of_linear_expression_lb > current_min_ub) {
min_of_linear_expression_lb = current_min_ub + 1;
}
if (min_of_linear_expression_lb > integer_trail_->LowerBound(min_var_)) {
integer_reason_.clear();
for (int i = 0; i < exprs_.size(); ++i) {
const IntegerValue slack = expr_lbs_[i] - min_of_linear_expression_lb;
integer_trail_->AppendRelaxedLinearReason(
slack, exprs_[i].coeffs, exprs_[i].vars, &integer_reason_);
}
if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(
min_var_, min_of_linear_expression_lb),
{}, integer_reason_)) {
return false;
}
}
// Propagation b) ub(min) >= ub(MIN(vars)) and we can't propagate anything
// here unless there is just one possible expression 'e' that can be the min:
// for all u != e, lb(u) > ub(min);
// In this case, ub(min) >= ub(e).
if (num_intervals_that_can_be_min == 1) {
const IntegerValue ub_of_only_candidate =
LinExprUpperBound(exprs_[last_possible_min_interval], *integer_trail_);
if (current_min_ub < ub_of_only_candidate) {
integer_reason_.clear();
// The reason is that all the other interval start after current_min_ub.
// And that min_ub has its current value.
integer_reason_.push_back(min_ub_literal);
for (int i = 0; i < exprs_.size(); ++i) {
if (i == last_possible_min_interval) continue;
const IntegerValue slack = expr_lbs_[i] - (current_min_ub + 1);
integer_trail_->AppendRelaxedLinearReason(
slack, exprs_[i].coeffs, exprs_[i].vars, &integer_reason_);
}
std::vector<Literal> literal_reason;
integer_trail_->MergeReasonInto(integer_reason_, &literal_reason);
auto constraint = absl::make_unique<IntegerSumLE>(
literal_reason, exprs_[last_possible_min_interval].vars,
exprs_[last_possible_min_interval].coeffs,
current_min_ub - exprs_[last_possible_min_interval].offset, model_);
constraint->Propagate();
ub_constraints_.emplace_back(std::move(constraint));
rev_ub_constraints_size_ = ub_constraints_.size();
}
}
return true;
}
void LinMinPropagator::RegisterWith(GenericLiteralWatcher* watcher) {
const int id = watcher->Register(this);
for (const LinearExpression& expr : exprs_) {
for (int i = 0; i < expr.vars.size(); ++i) {
const IntegerVariable& var = expr.vars[i];
const IntegerValue coeff = expr.coeffs[i];
if (coeff > 0) {
watcher->WatchLowerBound(var, id);
} else {
watcher->WatchUpperBound(var, id);
}
}
}
watcher->RegisterReversibleInt(id, &rev_ub_constraints_size_);
watcher->WatchUpperBound(min_var_, id);
}
PositiveProductPropagator::PositiveProductPropagator(
IntegerVariable a, IntegerVariable b, IntegerVariable p,
IntegerTrail* integer_trail)

View File

@@ -26,6 +26,7 @@
#include "ortools/sat/precedences.h"
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_solver.h"
#include "ortools/util/rev.h"
namespace operations_research {
namespace sat {
@@ -142,6 +143,57 @@ class MinPropagator : public PropagatorInterface {
DISALLOW_COPY_AND_ASSIGN(MinPropagator);
};
// Helper struct to model linear expression for lin_min/lin_max constraints. The
// canonical expression should only contain positive coefficients.
struct LinearExpression {
std::vector<IntegerVariable> vars;
std::vector<IntegerValue> coeffs;
IntegerValue offset = IntegerValue(0);
};
// Returns the same expression in the canonical form (all positive
// coefficients).
LinearExpression CanonicalizeExpr(const LinearExpression& expr);
// Returns lower bound of linear expression using variable bounds of the
// variables in expression. Assumes Canonical expression (all positive
// coefficients).
IntegerValue LinExprLowerBound(const LinearExpression& expr,
const IntegerTrail& integer_trail);
// Returns upper bound of linear expression using variable bounds of the
// variables in expression. Assumes Canonical expression (all positive
// coefficients).
IntegerValue LinExprUpperBound(const LinearExpression& expr,
const IntegerTrail& integer_trail);
// Same as MinPropagator except this works on min = MIN(exprs) where exprs are
// linear expressions. It uses IntegerSumLE to propagate bounds on the exprs.
// Assumes Canonical expressions (all positive coefficients).
// TODO(user): Remove the dependency on IntegerSumLE.
class LinMinPropagator : public PropagatorInterface {
public:
LinMinPropagator(const std::vector<LinearExpression>& exprs,
IntegerVariable min_var, Model* model);
LinMinPropagator(const LinMinPropagator&) = delete;
LinMinPropagator& operator=(const LinMinPropagator&) = delete;
bool Propagate() final;
void RegisterWith(GenericLiteralWatcher* watcher);
private:
const std::vector<LinearExpression> exprs_;
const IntegerVariable min_var_;
std::vector<IntegerValue> expr_lbs_;
Model* model_;
IntegerTrail* integer_trail_;
int rev_ub_constraints_size_ = 0;
std::vector<std::unique_ptr<IntegerSumLE>> ub_constraints_;
std::vector<IntegerLiteral> integer_reason_;
};
// Propagates a * b = c. Basic version, we don't extract any special cases, and
// we only propagates the bounds.
//
@@ -557,6 +609,59 @@ inline std::function<void(Model*)> IsEqualToMinOf(
};
}
// Expresses the fact that an existing integer variable is equal to the minimum
// of linear expressions. Assumes Canonical expressions (all positive
// coefficients).
inline std::function<void(Model*)> IsEqualToMinOf(
const LinearExpression& min_expr,
const std::vector<LinearExpression>& exprs) {
return [=](Model* model) {
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
IntegerVariable min_var;
if (min_expr.vars.size() == 1 &&
std::abs(min_expr.coeffs[0].value()) == 1) {
if (min_expr.coeffs[0].value() == 1) {
min_var = min_expr.vars[0];
} else {
min_var = NegationOf(min_expr.vars[0]);
}
} else {
// Create a new variable if the expression is not just a single variable.
IntegerValue min_lb = LinExprLowerBound(min_expr, *integer_trail);
IntegerValue min_ub = LinExprUpperBound(min_expr, *integer_trail);
min_var = integer_trail->AddIntegerVariable(min_lb, min_ub);
// min_var = min_expr
std::vector<IntegerVariable> min_sum_vars = min_expr.vars;
std::vector<int64> min_sum_coeffs;
for (IntegerValue coeff : min_expr.coeffs) {
min_sum_coeffs.push_back(coeff.value());
}
min_sum_vars.push_back(min_var);
min_sum_coeffs.push_back(-1);
model->Add(FixedWeightedSum(min_sum_vars, min_sum_coeffs,
-min_expr.offset.value()));
}
for (const LinearExpression& expr : exprs) {
// min_var <= expr
std::vector<IntegerVariable> vars = expr.vars;
std::vector<int64> coeffs;
for (IntegerValue coeff : expr.coeffs) {
coeffs.push_back(coeff.value());
}
vars.push_back(min_var);
coeffs.push_back(-1);
model->Add(WeightedSumGreaterOrEqual(vars, coeffs, -expr.offset.value()));
}
LinMinPropagator* constraint = new LinMinPropagator(exprs, min_var, model);
constraint->RegisterWith(model->GetOrCreate<GenericLiteralWatcher>());
model->TakeOwnership(constraint);
};
}
// Expresses the fact that an existing integer variable is equal to the maximum
// of other integer variables.
inline std::function<void(Model*)> IsEqualToMaxOf(
@@ -653,7 +758,7 @@ inline std::function<void(Model*)> ProductConstraint(IntegerVariable a,
};
}
// Adds the constraint: a / b = d.
// Adds the constraint: a / b = c.
inline std::function<void(Model*)> DivisionConstraint(IntegerVariable a,
IntegerVariable b,
IntegerVariable c) {
@@ -666,7 +771,7 @@ inline std::function<void(Model*)> DivisionConstraint(IntegerVariable a,
};
}
// Adds the constraint: a / b = d where b is a constant.
// Adds the constraint: a / b = c where b is a constant.
inline std::function<void(Model*)> FixedDivisionConstraint(IntegerVariable a,
IntegerValue b,
IntegerVariable c) {

View File

@@ -213,5 +213,32 @@ void MakeAllVariablesPositive(LinearConstraint* constraint) {
}
}
// TODO(user): it would be better if LinearConstraint natively supported
// term and not two separated vectors. Fix?
//
// TODO(user): This is really similar to CleanTermsAndFillConstraint(), maybe
// we should just make the later switch negative variable to positive ones to
// avoid an extra linear scan on each new cuts.
void CanonicalizeConstraint(LinearConstraint* ct) {
std::vector<std::pair<IntegerVariable, IntegerValue>> terms;
const int size = ct->vars.size();
for (int i = 0; i < size; ++i) {
if (VariableIsPositive(ct->vars[i])) {
terms.push_back({ct->vars[i], ct->coeffs[i]});
} else {
terms.push_back({NegationOf(ct->vars[i]), -ct->coeffs[i]});
}
}
std::sort(terms.begin(), terms.end());
ct->vars.clear();
ct->coeffs.clear();
for (const auto& term : terms) {
ct->vars.push_back(term.first);
ct->coeffs.push_back(term.second);
}
}
} // namespace sat
} // namespace operations_research

View File

@@ -148,6 +148,12 @@ void CleanTermsAndFillConstraint(
std::vector<std::pair<IntegerVariable, IntegerValue>>* terms,
LinearConstraint* constraint);
// Sorts the terms and makes all IntegerVariable positive. This assumes that a
// variable or its negation only appear once.
//
// Note that currently this allocates some temporary memory.
void CanonicalizeConstraint(LinearConstraint* ct);
} // namespace sat
} // namespace operations_research

View File

@@ -38,29 +38,6 @@ size_t ComputeHashOfTerms(const LinearConstraint& ct) {
return hash;
}
// TODO(user): it would be better if LinearConstraint natively supported
// term and not two separated vectors. Fix?
void CanonicalizeConstraint(LinearConstraint* ct) {
std::vector<std::pair<IntegerVariable, IntegerValue>> terms;
const int size = ct->vars.size();
for (int i = 0; i < size; ++i) {
if (VariableIsPositive(ct->vars[i])) {
terms.push_back({ct->vars[i], ct->coeffs[i]});
} else {
terms.push_back({NegationOf(ct->vars[i]), -ct->coeffs[i]});
}
}
std::sort(terms.begin(), terms.end());
ct->vars.clear();
ct->coeffs.clear();
for (const auto& term : terms) {
ct->vars.push_back(term.first);
ct->coeffs.push_back(term.second);
}
}
} // namespace
LinearConstraintManager::~LinearConstraintManager() {
@@ -132,7 +109,7 @@ bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
// to detect duplicate constraints and merge bounds. This is also relevant if
// we regenerate identical cuts for some reason.
LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
LinearConstraint ct) {
LinearConstraint ct, bool* added) {
CHECK(!ct.vars.empty());
SimplifyConstraint(&ct);
DivideByGCD(&ct);
@@ -145,25 +122,28 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
const ConstraintIndex ct_index = equiv_constraints_[key];
if (constraint_infos_[ct_index].constraint.vars == ct.vars &&
constraint_infos_[ct_index].constraint.coeffs == ct.coeffs) {
if (added != nullptr) *added = false;
if (ct.lb > constraint_infos_[ct_index].constraint.lb) {
if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
constraint_infos_[ct_index].constraint.lb = ct.lb;
if (added != nullptr) *added = true;
}
if (ct.ub < constraint_infos_[ct_index].constraint.ub) {
if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
constraint_infos_[ct_index].constraint.ub = ct.ub;
if (added != nullptr) *added = true;
}
++num_merged_constraints_;
return ct_index;
}
}
if (added != nullptr) *added = true;
const ConstraintIndex ct_index(constraint_infos_.size());
ConstraintInfo ct_info;
ct_info.constraint = std::move(ct);
ct_info.l2_norm = ComputeL2Norm(ct_info.constraint);
ct_info.is_in_lp = false;
ct_info.is_cut = false;
ct_info.objective_parallelism_computed = false;
ct_info.objective_parallelism = 0.0;
ct_info.inactive_count = 0;
@@ -227,19 +207,20 @@ void LinearConstraintManager::AddCut(
// Only add cut with sufficient efficacy.
if (violation / l2_norm < 1e-5) return;
VLOG(1) << "Cut '" << type_name << "'"
<< " size=" << ct.vars.size()
<< " max_magnitude=" << ComputeInfinityNorm(ct) << " norm=" << l2_norm
<< " violation=" << violation << " eff=" << violation / l2_norm;
// Add the constraint. We only mark the constraint as a cut if it is not an
// update of an already existing one.
const int64 prev_size = constraint_infos_.size();
const ConstraintIndex ct_index = Add(std::move(ct));
if (prev_size + 1 == constraint_infos_.size()) {
bool added = false;
const ConstraintIndex ct_index = Add(std::move(ct), &added);
if (added) {
VLOG(1) << "Cut '" << type_name << "'"
<< " size=" << constraint_infos_[ct_index].constraint.vars.size()
<< " max_magnitude="
<< ComputeInfinityNorm(constraint_infos_[ct_index].constraint)
<< " norm=" << l2_norm << " violation=" << violation
<< " eff=" << violation / l2_norm;
num_cuts_++;
type_to_num_cuts_[type_name]++;
constraint_infos_[ct_index].is_cut = true;
}
}

View File

@@ -45,7 +45,6 @@ class LinearConstraintManager {
struct ConstraintInfo {
LinearConstraint constraint;
double l2_norm;
bool is_cut;
int64 inactive_count;
double objective_parallelism;
bool objective_parallelism_computed;
@@ -66,9 +65,11 @@ class LinearConstraintManager {
// Add a new constraint to the manager. Note that we canonicalize constraints
// and merge the bounds of constraints with the same terms. We also perform
// basic preprocessing.
// basic preprocessing. If added is given, it will be set to true if this
// constraint was actually a new one and to false if it was dominated by an
// already existing one.
DEFINE_INT_TYPE(ConstraintIndex, int32);
ConstraintIndex Add(LinearConstraint ct);
ConstraintIndex Add(LinearConstraint ct, bool* added = nullptr);
// Same as Add(), but logs some information about the newly added constraint.
// Cuts are also handled slightly differently than normal constraints.

View File

@@ -1574,8 +1574,13 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint(
// If we add the row to the dense_terms, diff will indicate by how much
// |upper_bound - ImpliedLB(dense_terms)| will change. That correspond to
// increasing the multiplier by 1.
IntegerValue positive_diff = row_bound;
IntegerValue negative_diff = row_bound;
//
// At this stage, we are not sure computing sum coeff * bound will not
// overflow, so we use floating point numbers. It is fine to do so since
// this is not directly involved in the actual exact constraint generation:
// these variables are just used in an heuristic.
double positive_diff = ToDouble(row_bound);
double negative_diff = ToDouble(row_bound);
// TODO(user): we could relax a bit some of the condition and allow a sign
// change. It is just trickier to compute the diff when we allow such
@@ -1600,11 +1605,11 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint(
positive_limit = std::min(positive_limit, overflow_limit);
negative_limit = std::min(negative_limit, overflow_limit);
if (coeff > 0) {
positive_diff -= coeff * lb;
negative_diff -= coeff * ub;
positive_diff -= ToDouble(coeff) * ToDouble(lb);
negative_diff -= ToDouble(coeff) * ToDouble(ub);
} else {
positive_diff -= coeff * ub;
negative_diff -= coeff * lb;
positive_diff -= ToDouble(coeff) * ToDouble(ub);
negative_diff -= ToDouble(coeff) * ToDouble(lb);
}
continue;
}
@@ -1629,20 +1634,24 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint(
// This is how diff change.
const IntegerValue implied = current > 0 ? lb : ub;
if (implied != 0) {
positive_diff -= coeff * implied;
negative_diff -= coeff * implied;
positive_diff -= ToDouble(coeff) * ToDouble(implied);
negative_diff -= ToDouble(coeff) * ToDouble(implied);
}
}
// Only add a multiple of this row if it tighten the final constraint.
// The positive_diff/negative_diff are supposed to be integer modulo the
// double precision, so we only add a multiple if they seems far away from
// zero.
IntegerValue to_add(0);
if (positive_diff < 0 && positive_limit > 0) {
if (positive_diff <= -1.0 && positive_limit > 0) {
to_add = positive_limit;
}
if (negative_diff > 0 && negative_limit > 0) {
if (negative_diff >= 1.0 && negative_limit > 0) {
// Pick this if it is better than the positive sign.
if (to_add == 0 || IntTypeAbs(negative_limit * negative_diff) >
IntTypeAbs(positive_limit * positive_diff)) {
if (to_add == 0 ||
std::abs(ToDouble(negative_limit) * negative_diff) >
std::abs(ToDouble(positive_limit) * positive_diff)) {
to_add = -negative_limit;
}
}

View File

@@ -61,7 +61,7 @@ class Logger {
explicit Logger(LogBehavior v) : use_stdout_(v == STDOUT_LOG) {}
void Log(const std::string& message) {
if (use_stdout_) {
printf("%s\n", message.c_str());
absl::PrintF("%s\n", message);
} else {
LOG(INFO) << message;
}

View File

@@ -50,7 +50,8 @@ void PresolveContext::AddImplication(int a, int b) {
}
// b => x in [lb, ub].
void PresolveContext::AddImplyInDomain(int b, int x, const Domain& domain) {
ConstraintProto* PresolveContext::AddImplyInDomain(int b, int x,
const Domain& domain) {
ConstraintProto* const imply = working_model->add_constraints();
// Doing it like this seems to use slightly less memory.
@@ -60,6 +61,7 @@ void PresolveContext::AddImplyInDomain(int b, int x, const Domain& domain) {
mutable_linear->mutable_vars()->Resize(1, x);
mutable_linear->mutable_coeffs()->Resize(1, 1);
FillDomainInProto(domain, mutable_linear);
return imply;
}
bool PresolveContext::DomainIsEmpty(int ref) const {
@@ -424,9 +426,16 @@ void PresolveContext::InsertVarValueEncoding(int literal, int ref,
}
}
} else {
AddImplyInDomain(literal, var, Domain(var_value));
AddImplyInDomain(NegatedRef(literal), var,
Domain(var_value).Complement());
VLOG(2) << "Insert lit(" << literal << ") <=> var(" << var
<< ") == " << value;
const std::pair<int, int64> key{var, var_value};
eq_half_encoding[key][literal] =
AddImplyInDomain(literal, var, Domain(var_value));
VLOG(2) << " imply_eq = " << eq_half_encoding[key][literal];
neq_half_encoding[key][NegatedRef(literal)] = AddImplyInDomain(
NegatedRef(literal), var, Domain(var_value).Complement());
VLOG(2) << " imply_neq = "
<< neq_half_encoding[key][NegatedRef(literal)];
}
} else {
const int previous_literal = insert.first->second;
@@ -450,40 +459,38 @@ bool PresolveContext::InsertHalfVarValueEncoding(int literal, int var,
const auto& new_info = direct_map.insert(std::make_pair(literal, ct));
if (new_info.second) {
VLOG(2) << "Collect lit(" << literal << ") implies var(" << var
<< (imply_eq ? ") == " : ") != ") << value;
<< (imply_eq ? ") == " : ") != ") << value << " [" << ct << "]";
UpdateRuleStats("variables: detect half reified value encoding");
if (gtl::ContainsKey(other_map, NegatedRef(literal))) {
const int ref_lit = imply_eq ? literal : NegatedRef(literal);
const auto insert_encoding_status =
encoding.insert(std::make_pair(std::make_pair(var, value), ref_lit));
if (insert_encoding_status.second) {
VLOG(2) << "Detect and store lit(" << ref_lit << ") <=> var(" << var
<< ") == " << value;
UpdateRuleStats("variables: detect fully reified value encoding");
} else if (ref_lit != insert_encoding_status.first->second) {
const int previous_lit = insert_encoding_status.first->second;
VLOG(2) << "Detect duplicate encoding lit(" << ref_lit << ") == lit("
<< previous_lit << ") <=> var(" << var << ") == " << value;
StoreBooleanEqualityRelation(previous_lit, ref_lit);
const int imply_eq_literal = imply_eq ? literal : NegatedRef(literal);
const AffineRelation::Relation r = GetAffineRelation(previous_lit);
if (r.representative != previous_lit) {
const int new_ref_lit =
r.coeff > 0 ? r.representative : NegatedRef(r.representative);
VLOG(2)
<< "Different representative of previous literal, updating to "
<< new_ref_lit;
const auto insert_encoding_status =
encoding.insert(std::make_pair(key, imply_eq_literal));
if (insert_encoding_status.second) {
VLOG(2) << "Detect and store lit(" << imply_eq_literal << ") <=> var("
<< var << ") == " << value;
UpdateRuleStats("variables: detect fully reified value encoding");
encoding[key] = imply_eq_literal;
} else if (imply_eq_literal != insert_encoding_status.first->second) {
const int previous_imply_eq_literal =
insert_encoding_status.first->second;
VLOG(2) << "Detect duplicate encoding lit(" << imply_eq_literal
<< ") == lit(" << previous_imply_eq_literal << ") <=> var("
<< var << ") == " << value;
StoreBooleanEqualityRelation(imply_eq_literal,
previous_imply_eq_literal);
// Update reference literal.
const AffineRelation::Relation r =
GetAffineRelation(previous_imply_eq_literal);
const int new_ref_lit =
r.coeff > 0 ? r.representative : NegatedRef(r.representative);
if (new_ref_lit != previous_imply_eq_literal) {
VLOG(2) << "Updating reference encoding literal to " << new_ref_lit;
encoding[key] = new_ref_lit;
}
// Remove the two half encoding constraints from the model.
VLOG(2) << "Delete " << ct->DebugString();
ct->Clear();
direct_map.erase(literal);
VLOG(2) << "Delete " << other_map[NegatedRef(literal)]->DebugString();
other_map[NegatedRef(literal)]->Clear();
other_map.erase(NegatedRef(literal));
UpdateRuleStats(
"variables: merge equivalent var value encoding literals");
}
@@ -534,6 +541,8 @@ int PresolveContext::GetOrCreateVarValueEncoding(int ref, int64 value) {
const int64 var_min = MinOf(var);
const int64 var_max = MaxOf(var);
if (domains[var].Size() == 2) {
VLOG(2) << "GetOrCreate var(" << var << ") {" << var_min << ", " << var_max
<< "}";
// Checks if the other value is already encoded.
const int64 other_value = var_value == var_min ? var_max : var_min;
const std::pair<int, int64> other_key{var, other_value};
@@ -724,7 +733,8 @@ void PresolveContext::SubstituteVariableInObjective(
// We can only "easily" substitute if the objective coefficient is a multiple
// of the one in the constraint.
const int64 coeff_in_objective = gtl::FindOrDie(objective_map, var_in_equality);
const int64 coeff_in_objective =
gtl::FindOrDie(objective_map, var_in_equality);
CHECK_NE(coeff_in_equality, 0);
CHECK_EQ(coeff_in_objective % coeff_in_equality, 0);
const int64 multiplier = coeff_in_objective / coeff_in_equality;
@@ -765,6 +775,11 @@ void PresolveContext::SubstituteVariableInObjective(
// Tricky: The objective domain is without the offset, so we need to shift it
objective_offset += static_cast<double>(offset.Min());
objective_domain = objective_domain.AdditionWith(Domain(-offset.Min()));
// Because we can assume that the constraint we used was constraining
// (otherwise it would have been removed), the objective domain should be now
// constraining.
objective_domain_is_constraining = true;
}
void PresolveContext::WriteObjectiveToProto() {

View File

@@ -48,7 +48,7 @@ struct PresolveContext {
void AddImplication(int a, int b);
// b => x in [lb, ub].
void AddImplyInDomain(int b, int x, const Domain& domain);
ConstraintProto* AddImplyInDomain(int b, int x, const Domain& domain);
// Helpers to query the current domain of a variable.
bool DomainIsEmpty(int ref) const;