works on CP-SAT and cuts

This commit is contained in:
Laurent Perron
2019-07-18 11:36:47 -07:00
parent d216e3caab
commit f32e012aed
10 changed files with 156 additions and 123 deletions

View File

@@ -57,6 +57,11 @@ message IntegerVariableProto {
// - max == domain.back();
// - for all i < n : min_i <= max_i
// - for all i < n-1 : max_i + 1 < min_{i+1}.
//
// Note that we check at validation that a variable domain is small enough so
// that we don't run into integer overflow in our algorithms. Because of that,
// you cannot just have "unbounded" variable like [0, kint64max] and should
// try to specify tighter domains.
repeated int64 domain = 2;
}

View File

@@ -81,16 +81,26 @@ std::string ValidateIntegerVariable(const CpModelProto& model, int v) {
ProtobufShortDebugString(proto));
}
// We do compute ub - lb in some place in the code and do not want to deal
// with overflow everywhere. This seems like a reasonable precondition anyway.
// Internally, we often take the negation of a domain, and we also want to
// have sentinel values greater than the min/max of a variable domain, so
// the domain must fall in [kint64min + 2, kint64max - 1].
const int64 lb = proto.domain(0);
const int64 ub = proto.domain(proto.domain_size() - 1);
if (lb < kint64min + 2 || ub > kint64max - 1) {
return absl::StrCat(
"var #", v, " domain do not fall in [kint64min + 2, kint64max - 1]. ",
ProtobufShortDebugString(proto));
}
// We do compute ub - lb in some place in the code and do not want to deal
// with overflow everywhere. This seems like a reasonable precondition anyway.
if (lb < 0 && lb + kint64max < ub) {
return absl::StrCat(
"var #", v,
" has a domain that is too large, i.e. |UB - LB| overflow an int64: ",
ProtobufShortDebugString(proto));
}
return "";
}

View File

@@ -28,7 +28,10 @@ namespace sat {
// fails at the first error and returns a human-readable description of the
// issue.
//
// TODO(user): Add any needed overflow validation.
// TODO(user): Add any needed overflow validation because we are far from
// exhaustive. We could also run a small presolve that tighten variable bounds
// before the overflow check to facilitate the lives of our users, but it is a
// some work to put in place.
std::string ValidateCpModel(const CpModelProto& model);
// Verifies that the given variable assignment is a feasible solution of the

View File

@@ -430,43 +430,33 @@ void TryToAddCutGenerators(const CpModelProto& model_proto,
std::vector<Literal> literals = mapping->Literals(ct.circuit().literals());
const int num_nodes = ReindexArcs(&tails, &heads, &literals);
std::vector<IntegerVariable> vars;
vars.reserve(literals.size());
for (const Literal& literal : literals) {
vars.push_back(m->Add(NewIntegerVariableFromLiteral(literal)));
}
relaxation->cut_generators.push_back(
CreateStronglyConnectedGraphCutGenerator(num_nodes, tails, heads,
vars));
literals, m));
}
if (ct.constraint_case() == ConstraintProto::ConstraintCase::kRoutes &&
linearization_level > 1) {
std::vector<int> tails;
std::vector<int> heads;
std::vector<IntegerVariable> vars;
int num_nodes = 0;
auto* encoder = m->GetOrCreate<IntegerEncoder>();
for (int i = 0; i < ct.routes().tails_size(); ++i) {
const IntegerVariable var =
encoder->GetLiteralView(mapping->Literal(ct.routes().literals(i)));
if (var == kNoIntegerVariable) return;
std::vector<int> tails(ct.routes().tails().begin(),
ct.routes().tails().end());
std::vector<int> heads(ct.routes().heads().begin(),
ct.routes().heads().end());
std::vector<Literal> literals = mapping->Literals(ct.routes().literals());
vars.push_back(var);
tails.push_back(ct.routes().tails(i));
heads.push_back(ct.routes().heads(i));
int num_nodes = 0;
for (int i = 0; i < ct.routes().tails_size(); ++i) {
num_nodes = std::max(num_nodes, 1 + ct.routes().tails(i));
num_nodes = std::max(num_nodes, 1 + ct.routes().heads(i));
}
if (ct.routes().demands().empty() || ct.routes().capacity() == 0) {
relaxation->cut_generators.push_back(
CreateStronglyConnectedGraphCutGenerator(num_nodes, tails, heads,
vars));
literals, m));
} else {
const std::vector<int64> demands(ct.routes().demands().begin(),
ct.routes().demands().end());
relaxation->cut_generators.push_back(CreateCVRPCutGenerator(
num_nodes, tails, heads, vars, demands, ct.routes().capacity()));
relaxation->cut_generators.push_back(
CreateCVRPCutGenerator(num_nodes, tails, heads, literals, demands,
ct.routes().capacity(), m));
}
}
if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntProd) {
@@ -529,11 +519,16 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto,
auto* mapping = m->GetOrCreate<CpModelMapping>();
auto* encoder = m->GetOrCreate<IntegerEncoder>();
auto* trail = m->GetOrCreate<Trail>();
for (const auto& ct : model_proto.constraints()) {
// Make sure the literal from a circuit constraint always have a view.
// Make sure the literals from a circuit constraint always have a view.
if (ct.constraint_case() == ConstraintProto::ConstraintCase::kCircuit) {
for (const int ref : ct.circuit().literals()) {
m->Add(NewIntegerVariableFromLiteral(mapping->Literal(ref)));
const Literal l = mapping->Literal(ref);
if (encoder->GetLiteralView(l) == kNoIntegerVariable &&
encoder->GetLiteralView(l.Negated()) == kNoIntegerVariable) {
m->Add(NewIntegerVariableFromLiteral(l));
}
}
}
@@ -546,8 +541,12 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto,
bool ok = true;
for (const int literal_ref : refs.literals) {
const Literal literal = mapping->Literal(literal_ref);
if (encoder->GetLiteralView(literal) == kNoIntegerVariable &&
encoder->GetLiteralView(literal.Negated()) == kNoIntegerVariable) {
if (trail->Assignment().LiteralIsAssigned(literal)) {
// Create a view to the constant 0 or 1.
m->Add(NewIntegerVariableFromLiteral(literal));
} else if (encoder->GetLiteralView(literal) == kNoIntegerVariable &&
encoder->GetLiteralView(literal.Negated()) ==
kNoIntegerVariable) {
ok = false;
break;
}
@@ -556,7 +555,6 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto,
TryToLinearizeConstraint(model_proto, ct, m, linearization_level,
&relaxation);
TryToAddCutGenerators(model_proto, ct, m, &relaxation);
}

View File

@@ -31,15 +31,6 @@ void LinearConstraintBuilder::AddTerm(IntegerVariable var, IntegerValue coeff) {
ABSL_MUST_USE_RESULT bool LinearConstraintBuilder::AddLiteralTerm(
Literal lit, IntegerValue coeff) {
if (assignment_.LiteralIsTrue(lit)) {
if (lb_ > kMinIntegerValue) lb_ -= coeff;
if (ub_ < kMaxIntegerValue) ub_ -= coeff;
return true;
}
if (assignment_.LiteralIsFalse(lit)) {
return true;
}
bool has_direct_view = encoder_.GetLiteralView(lit) != kNoIntegerVariable;
bool has_opposite_view =
encoder_.GetLiteralView(lit.Negated()) != kNoIntegerVariable;

View File

@@ -77,16 +77,14 @@ struct LinearConstraint {
};
// Allow to build a LinearConstraint while making sure there is no duplicate
// variables.
// variables. Note that we do not simplify literal/variable that are currently
// fixed here.
class LinearConstraintBuilder {
public:
// We support "sticky" kMinIntegerValue for lb and kMaxIntegerValue for ub
// for one-sided constraints.
LinearConstraintBuilder(const Model* model, IntegerValue lb, IntegerValue ub)
: assignment_(model->Get<Trail>()->Assignment()),
encoder_(*model->Get<IntegerEncoder>()),
lb_(lb),
ub_(ub) {}
: encoder_(*model->Get<IntegerEncoder>()), lb_(lb), ub_(ub) {}
// Adds var * coeff to the constraint.
void AddTerm(IntegerVariable var, IntegerValue coeff);
@@ -105,7 +103,6 @@ class LinearConstraintBuilder {
LinearConstraint Build();
private:
const VariablesAssignment& assignment_;
const IntegerEncoder& encoder_;
IntegerValue lb_;
IntegerValue ub_;

View File

@@ -800,7 +800,7 @@ bool LinearProgrammingConstraint::Propagate() {
const IntegerValue propagated_lb =
integer_trail_->LowerBound(objective_cp_);
if (approximate_new_lb > propagated_lb) {
VLOG(1) << "LP objective [ " << ToDouble(propagated_lb) << ", "
VLOG(2) << "LP objective [ " << ToDouble(propagated_lb) << ", "
<< ToDouble(integer_trail_->UpperBound(objective_cp_))
<< " ] approx_lb += "
<< ToDouble(approximate_new_lb - propagated_lb);
@@ -1467,22 +1467,20 @@ void AddOutgoingCut(int num_nodes, int subset_size,
const std::vector<bool>& in_subset,
const std::vector<int>& tails,
const std::vector<int>& heads,
const std::vector<IntegerVariable>& vars,
const std::vector<double>& var_lp_values,
const std::vector<Literal>& literals,
const std::vector<double>& literal_lp_values,
int64 rhs_lower_bound,
const gtl::ITIVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
LinearConstraint outgoing;
LinearConstraintManager* manager, Model* model) {
LinearConstraintBuilder outgoing(model, IntegerValue(rhs_lower_bound),
kMaxIntegerValue);
double sum_outgoing = 0.0;
outgoing.lb = IntegerValue(rhs_lower_bound);
outgoing.ub = kMaxIntegerValue;
// Add outgoing arcs, compute outgoing flow.
for (int i = 0; i < tails.size(); ++i) {
if (in_subset[tails[i]] && !in_subset[heads[i]]) {
sum_outgoing += var_lp_values[i];
outgoing.vars.push_back(vars[i]);
outgoing.coeffs.push_back(IntegerValue(1));
sum_outgoing += literal_lp_values[i];
CHECK(outgoing.AddLiteralTerm(literals[i], IntegerValue(1)));
}
}
@@ -1501,13 +1499,13 @@ void AddOutgoingCut(int num_nodes, int subset_size,
if (in_subset[tails[i]]) {
num_optional_nodes_in++;
if (optional_loop_in == -1 ||
var_lp_values[i] < var_lp_values[optional_loop_in]) {
literal_lp_values[i] < literal_lp_values[optional_loop_in]) {
optional_loop_in = i;
}
} else {
num_optional_nodes_out++;
if (optional_loop_out == -1 ||
var_lp_values[i] < var_lp_values[optional_loop_out]) {
literal_lp_values[i] < literal_lp_values[optional_loop_out]) {
optional_loop_out = i;
}
}
@@ -1517,32 +1515,32 @@ void AddOutgoingCut(int num_nodes, int subset_size,
// When all optionals of one side are excluded in lp solution, no cut.
if (num_optional_nodes_in == subset_size &&
(optional_loop_in == -1 ||
var_lp_values[optional_loop_in] > 1.0 - 1e-6)) {
literal_lp_values[optional_loop_in] > 1.0 - 1e-6)) {
return;
}
if (num_optional_nodes_out == num_nodes - subset_size &&
(optional_loop_out == -1 ||
var_lp_values[optional_loop_out] > 1.0 - 1e-6)) {
literal_lp_values[optional_loop_out] > 1.0 - 1e-6)) {
return;
}
// There is no mandatory node in subset, add optional_loop_in.
if (num_optional_nodes_in == subset_size) {
outgoing.vars.push_back(vars[optional_loop_in]);
outgoing.coeffs.push_back(IntegerValue(1));
sum_outgoing += var_lp_values[optional_loop_in];
CHECK(
outgoing.AddLiteralTerm(literals[optional_loop_in], IntegerValue(1)));
sum_outgoing += literal_lp_values[optional_loop_in];
}
// There is no mandatory node out of subset, add optional_loop_out.
if (num_optional_nodes_out == num_nodes - subset_size) {
outgoing.vars.push_back(vars[optional_loop_out]);
outgoing.coeffs.push_back(IntegerValue(1));
sum_outgoing += var_lp_values[optional_loop_out];
CHECK(outgoing.AddLiteralTerm(literals[optional_loop_out],
IntegerValue(1)));
sum_outgoing += literal_lp_values[optional_loop_out];
}
}
if (sum_outgoing < rhs_lower_bound - 1e-6) {
manager->AddCut(outgoing, "Circuit", lp_values);
manager->AddCut(outgoing.Build(), "Circuit", lp_values);
}
}
@@ -1556,10 +1554,10 @@ void AddOutgoingCut(int num_nodes, int subset_size,
// the asymmetric case.
void SeparateSubtourInequalities(
int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
const std::vector<IntegerVariable>& vars,
const std::vector<Literal>& literals,
const gtl::ITIVector<IntegerVariable, double>& lp_values,
absl::Span<const int64> demands, int64 capacity,
LinearConstraintManager* manager) {
LinearConstraintManager* manager, Model* model) {
if (num_nodes <= 2) return;
// We will collect only the arcs with a positive lp_values to speed up some
@@ -1572,11 +1570,20 @@ void SeparateSubtourInequalities(
std::vector<Arc> relevant_arcs;
// Sort the arcs by non-increasing lp_values.
std::vector<double> literal_lp_values(literals.size());
std::vector<std::pair<double, int>> arc_by_decreasing_lp_values;
std::vector<double> var_lp_values;
for (int i = 0; i < vars.size(); ++i) {
const double lp_value = lp_values[vars[i]];
var_lp_values.push_back(lp_value);
auto* encoder = model->GetOrCreate<IntegerEncoder>();
for (int i = 0; i < literals.size(); ++i) {
double lp_value;
const IntegerVariable direct_view = encoder->GetLiteralView(literals[i]);
if (direct_view != kNoIntegerVariable) {
lp_value = lp_values[direct_view];
} else {
lp_value =
1.0 - lp_values[encoder->GetLiteralView(literals[i].Negated())];
}
literal_lp_values[i] = lp_value;
if (lp_value < 1e-6) continue;
relevant_arcs.push_back({tails[i], heads[i], lp_value});
arc_by_decreasing_lp_values.push_back({lp_value, i});
@@ -1751,9 +1758,10 @@ void SeparateSubtourInequalities(
// Add a cut if the current outgoing flow is not enough.
if (outgoing_flow < min_outgoing_flow - 1e-6) {
AddOutgoingCut(num_nodes, subset.size(), in_subset, tails, heads, vars,
var_lp_values,
/*rhs_lower_bound=*/min_outgoing_flow, lp_values, manager);
AddOutgoingCut(num_nodes, subset.size(), in_subset, tails, heads,
literals, literal_lp_values,
/*rhs_lower_bound=*/min_outgoing_flow, lp_values, manager,
model);
}
// Sparse clean up.
@@ -1761,20 +1769,42 @@ void SeparateSubtourInequalities(
}
}
namespace {
// Returns for each literal its integer view, or the view of its negation.
std::vector<IntegerVariable> GetAssociatedVariables(
const std::vector<Literal>& literals, Model* model) {
auto* encoder = model->GetOrCreate<IntegerEncoder>();
std::vector<IntegerVariable> result;
for (const Literal l : literals) {
const IntegerVariable direct_view = encoder->GetLiteralView(l);
if (direct_view != kNoIntegerVariable) {
result.push_back(direct_view);
} else {
result.push_back(encoder->GetLiteralView(l.Negated()));
DCHECK_NE(result.back(), kNoIntegerVariable);
}
}
return result;
}
} // namespace
// We use a basic algorithm to detect components that are not connected to the
// rest of the graph in the LP solution, and add cuts to force some arcs to
// enter and leave this component from outside.
CutGenerator CreateStronglyConnectedGraphCutGenerator(
int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
const std::vector<IntegerVariable>& vars) {
const std::vector<Literal>& literals, Model* model) {
CutGenerator result;
result.vars = vars;
result.vars = GetAssociatedVariables(literals, model);
result.generate_cuts =
[num_nodes, tails, heads, vars](
[num_nodes, tails, heads, literals, model](
const gtl::ITIVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
SeparateSubtourInequalities(num_nodes, tails, heads, vars, lp_values,
/*demands=*/{}, /*capacity=*/0, manager);
SeparateSubtourInequalities(
num_nodes, tails, heads, literals, lp_values,
/*demands=*/{}, /*capacity=*/0, manager, model);
};
return result;
}
@@ -1782,17 +1812,18 @@ CutGenerator CreateStronglyConnectedGraphCutGenerator(
CutGenerator CreateCVRPCutGenerator(int num_nodes,
const std::vector<int>& tails,
const std::vector<int>& heads,
const std::vector<IntegerVariable>& vars,
const std::vector<Literal>& literals,
const std::vector<int64>& demands,
int64 capacity) {
int64 capacity, Model* model) {
CutGenerator result;
result.vars = vars;
result.vars = GetAssociatedVariables(literals, model);
result.generate_cuts =
[num_nodes, tails, heads, demands, capacity, vars](
[num_nodes, tails, heads, demands, capacity, literals, model](
const gtl::ITIVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
SeparateSubtourInequalities(num_nodes, tails, heads, vars, lp_values,
demands, capacity, manager);
SeparateSubtourInequalities(num_nodes, tails, heads, literals,
lp_values, demands, capacity, manager,
model);
};
return result;
}

View File

@@ -428,7 +428,7 @@ class LinearProgrammingConstraintCollection
// we do not add any cuts for components of size 1.
CutGenerator CreateStronglyConnectedGraphCutGenerator(
int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
const std::vector<IntegerVariable>& vars);
const std::vector<Literal>& literals, Model* model);
// Almost the same as CreateStronglyConnectedGraphCutGenerator() but for each
// components, computes the demand needed to serves it, and depending on whether
@@ -437,9 +437,9 @@ CutGenerator CreateStronglyConnectedGraphCutGenerator(
CutGenerator CreateCVRPCutGenerator(int num_nodes,
const std::vector<int>& tails,
const std::vector<int>& heads,
const std::vector<IntegerVariable>& vars,
const std::vector<Literal>& literals,
const std::vector<int64>& demands,
int64 capacity);
int64 capacity, Model* model);
} // namespace sat
} // namespace operations_research

View File

@@ -176,8 +176,8 @@ class LinearExpr(object):
while to_process: # Flatten to avoid recursion.
expr, coef = to_process.pop()
if isinstance(expr, _ProductCst):
to_process.append((expr.Expression(),
coef * expr.Coefficient()))
to_process.append(
(expr.Expression(), coef * expr.Coefficient()))
elif isinstance(expr, _SumArray):
for e in expr.Expressions():
to_process.append((e, coef))
@@ -423,8 +423,8 @@ class _ScalProd(LinearExpr):
def __repr__(self):
return 'ScalProd([{}], [{}], {})'.format(
', '.join(map(repr, self.__expressions)), ', '.join(
map(repr, self.__coefficients)), self.__constant)
', '.join(map(repr, self.__expressions)),
', '.join(map(repr, self.__coefficients)), self.__constant)
def Expressions(self):
return self.__expressions
@@ -663,14 +663,14 @@ class IntervalVar(object):
if self.__ct.enforcement_literal:
return '%s(start = %s, size = %s, end = %s, is_present = %s)' % (
self.__ct.name, ShortName(self.__model, interval.start),
ShortName(self.__model, interval.size),
ShortName(self.__model, interval.end),
ShortName(self.__model,
interval.size), ShortName(self.__model, interval.end),
ShortName(self.__model, self.__ct.enforcement_literal[0]))
else:
return '%s(start = %s, size = %s, end = %s)' % (
self.__ct.name, ShortName(self.__model, interval.start),
ShortName(self.__model, interval.size),
ShortName(self.__model, interval.end))
ShortName(self.__model,
interval.size), ShortName(self.__model, interval.end))
def Name(self):
return self.__ct.name
@@ -1548,8 +1548,8 @@ def EvaluateBooleanExpression(literal, solution):
else:
return not solution.solution[-index - 1]
else:
raise TypeError(
'Cannot interpret %s as a boolean expression.' % literal)
raise TypeError('Cannot interpret %s as a boolean expression.' %
literal)
class CpSolver(object):
@@ -1717,8 +1717,8 @@ class CpSolverSolutionCallback(pywrapsat.SolutionCallback):
index = lit.Index()
return self.SolutionBooleanValue(index)
else:
raise TypeError(
'Cannot interpret %s as a boolean expression.' % lit)
raise TypeError('Cannot interpret %s as a boolean expression.' %
lit)
def Value(self, expression):
"""Evaluates an linear expression in the current solution.
@@ -1742,8 +1742,8 @@ class CpSolverSolutionCallback(pywrapsat.SolutionCallback):
while to_process:
expr, coef = to_process.pop()
if isinstance(expr, _ProductCst):
to_process.append((expr.Expression(),
coef * expr.Coefficient()))
to_process.append(
(expr.Expression(), coef * expr.Coefficient()))
elif isinstance(expr, _SumArray):
for e in expr.Expressions():
to_process.append((e, coef))
@@ -1755,8 +1755,8 @@ class CpSolverSolutionCallback(pywrapsat.SolutionCallback):
elif isinstance(expr, IntVar):
value += coef * self.SolutionIntegerValue(expr.Index())
elif isinstance(expr, _NotBooleanVariable):
value += coef * (
1 - self.SolutionIntegerValue(expr.Not().Index()))
value += coef * (1 -
self.SolutionIntegerValue(expr.Not().Index()))
return value

View File

@@ -71,11 +71,10 @@ def DisplayJobshop(starts, durations, machines, name):
for i in all_jobs:
for j in all_machines:
df.append(
dict(
Task='Resource%i' % machines[i][j],
Start=ToDate(starts[i][j]),
Finish=ToDate(starts[i][j] + durations[i][j]),
Resource='Job%i' % i))
dict(Task='Resource%i' % machines[i][j],
Start=ToDate(starts[i][j]),
Finish=ToDate(starts[i][j] + durations[i][j]),
Resource='Job%i' % i))
sorted_df = sorted(df, key=lambda k: k['Task'])
@@ -87,15 +86,14 @@ def DisplayJobshop(starts, durations, machines, name):
for i in all_jobs:
colors['Job%i' % i] = cm.RandomColor()
fig = ff.create_gantt(
sorted_df,
colors=colors,
index_col='Resource',
title=name,
show_colorbar=False,
showgrid_x=True,
showgrid_y=True,
group_tasks=True)
fig = ff.create_gantt(sorted_df,
colors=colors,
index_col='Resource',
title=name,
show_colorbar=False,
showgrid_x=True,
showgrid_y=True,
group_tasks=True)
pyo.iplot(fig)
@@ -144,8 +142,8 @@ class SvgWrapper(object):
self.__dwg.line((o, y), (self.__sizex * s + o, y), stroke='black'))
for i in range(0, int(self.__sizex) + 1, step):
self.__dwg.add(
self.__dwg.line(
(o + i * s, y - dy), (o + i * s, y + dy), stroke='black'))
self.__dwg.line((o + i * s, y - dy), (o + i * s, y + dy),
stroke='black'))
def AddYScale(self, step=1):
"""Add an scale on the y axis."""
@@ -157,8 +155,8 @@ class SvgWrapper(object):
self.__dwg.line((x, o), (x, self.__sizey * s + o), stroke='black'))
for i in range(0, int(self.__sizey) + 1, step):
self.__dwg.add(
self.__dwg.line(
(x - dx, i * s + o), (x + dx, i * s + o), stroke='black'))
self.__dwg.line((x - dx, i * s + o), (x + dx, i * s + o),
stroke='black'))
def AddTitle(self, title):
"""Add a title to the drawing."""