From 9d1803cee253b30e475bfcd7d441aa20be0fe93d Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Thu, 7 Apr 2022 17:59:22 +0200 Subject: [PATCH] [CP-SAT] bugfixes; more work on cuts on graphs --- ortools/sat/cp_model_presolve.cc | 82 ++++ ortools/sat/integer.h | 2 + ortools/sat/linear_programming_constraint.cc | 376 ++++++++++++++----- ortools/sat/linear_programming_constraint.h | 49 +++ ortools/sat/scheduling_constraints.cc | 10 +- 5 files changed, 423 insertions(+), 96 deletions(-) diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 972519fa86..0d287baaaa 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -1310,6 +1310,55 @@ bool CpModelPresolver::PresolveIntMod(ConstraintProto* ct) { const LinearExpressionProto expr = ct->int_mod().exprs(0); const LinearExpressionProto mod = ct->int_mod().exprs(1); + if (context_->MinOf(target) > 0) { + bool domain_changed = false; + if (!context_->IntersectDomainWith( + expr, Domain(0, std::numeric_limits::max()), + &domain_changed)) { + return false; + } + if (domain_changed) { + context_->UpdateRuleStats( + "int_mod: non negative target implies positive expression"); + } + } + + if (context_->MinOf(target) >= context_->MaxOf(mod) || + context_->MaxOf(target) <= -context_->MaxOf(mod)) { + return context_->NotifyThatModelIsUnsat( + "int_mod: incompatible target and mod"); + } + + if (context_->MaxOf(target) < 0) { + bool domain_changed = false; + if (!context_->IntersectDomainWith( + expr, Domain(std::numeric_limits::min(), 0), + &domain_changed)) { + return false; + } + if (domain_changed) { + context_->UpdateRuleStats( + "int_mod: non positive target implies negative expression"); + } + } + + if (context_->IsFixed(target) && context_->IsFixed(mod) && + context_->FixedValue(mod) > 1 && ct->enforcement_literal().empty() && + expr.vars().size() == 1) { + // We can intersect the domain of expr with {k * mod + target}. + const int64_t fixed_mod = context_->FixedValue(mod); + const int64_t fixed_target = context_->FixedValue(target); + + if (!context_->CanonicalizeAffineVariable(expr.vars(0), expr.coeffs(0), + fixed_mod, + fixed_target - expr.offset())) { + return false; + } + + context_->UpdateRuleStats("int_mod: fixed mod and target"); + return RemoveConstraint(ct); + } + bool domain_changed = false; if (!context_->IntersectDomainWith( target, @@ -8487,6 +8536,39 @@ CpSolverStatus CpModelPresolver::Presolve() { EncodeAllAffineRelations(); if (context_->ModelIsUnsat()) return InfeasibleStatus(); + // Final cleaning of the scheduling constraints. Probing could have set + // some literal to false, turning the intervals into empty constraints, while + // not having cleaned up the scheduling constraints. + { + const int num_constraints = context_->working_model->constraints_size(); + for (int c = 0; c < num_constraints; ++c) { + ConstraintProto* ct = context_->working_model->mutable_constraints(c); + switch (ct->constraint_case()) { + case ConstraintProto::kNoOverlap: { + // Filter out absent intervals. + if (PresolveNoOverlap(ct)) { + context_->UpdateConstraintVariableUsage(c); + } + break; + } + case ConstraintProto::kNoOverlap2D: + // Filter out absent intervals. + if (PresolveNoOverlap2D(c, ct)) { + context_->UpdateConstraintVariableUsage(c); + } + break; + case ConstraintProto::kCumulative: + // Filter out absent intervals. + if (PresolveCumulative(ct)) { + context_->UpdateConstraintVariableUsage(c); + } + break; + default: { + } + } + } + } + // The strategy variable indices will be remapped in ApplyVariableMapping() // but first we use the representative of the affine relations for the // variables that are not present anymore. diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index eb83ab6def..c61c29b3f7 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -280,6 +280,8 @@ struct AffineExpression { return ToDouble(coeff) * lp_values[var] + ToDouble(constant); } + bool IsConstant() const { return var == kNoIntegerVariable; } + const std::string DebugString() const { if (var == kNoIntegerVariable) return absl::StrCat(constant.value()); if (constant == 0) { diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 95c0c7534c..40f047e57c 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -2377,6 +2377,109 @@ void AddOutgoingCut( } // namespace +void GenerateInterestingSubsets(int num_nodes, + const std::vector>& arcs, + int min_subset_size, int stop_at_num_components, + std::vector* subset_data, + std::vector>* subsets) { + subset_data->resize(num_nodes); + subsets->clear(); + + // We will do a union-find by adding one by one the arc of the lp solution + // in the order above. Every intermediate set during this construction will + // be a candidate for a cut. + // + // In parallel to the union-find, to efficiently reconstruct these sets (at + // most num_nodes), we construct a "decomposition forest" of the different + // connected components. Note that we don't exploit any asymmetric nature of + // the graph here. This is exactly the algo 6.3 in the book above. + int num_components = num_nodes; + std::vector parent(num_nodes); + std::vector root(num_nodes); + for (int i = 0; i < num_nodes; ++i) { + parent[i] = i; + root[i] = i; + } + auto get_root_and_compress_path = [&root](int node) { + int r = node; + while (root[r] != r) r = root[r]; + while (root[node] != r) { + const int next = root[node]; + root[node] = r; + node = next; + } + return r; + }; + for (const auto& [initial_tail, initial_head] : arcs) { + if (num_components <= stop_at_num_components) break; + const int tail = get_root_and_compress_path(initial_tail); + const int head = get_root_and_compress_path(initial_head); + if (tail != head) { + // Update the decomposition forest, note that the number of nodes is + // growing. + const int new_node = parent.size(); + parent.push_back(new_node); + parent[head] = new_node; + parent[tail] = new_node; + --num_components; + + // It is important that the union-find representative is the same node. + root.push_back(new_node); + root[head] = new_node; + root[tail] = new_node; + } + } + + // For each node in the decomposition forest, try to add a cut for the set + // formed by the nodes and its children. To do that efficiently, we first + // order the nodes so that for each node in a tree, the set of children forms + // a consecutive span in the subset_data vector. This vector just lists the + // nodes in the "pre-order" graph traversal order. The Spans will point inside + // the subset_data vector, it is why we initialize it once and for all. + int new_size = 0; + { + std::vector> graph(parent.size()); + for (int i = 0; i < parent.size(); ++i) { + if (parent[i] != i) graph[parent[i]].push_back(i); + } + std::vector queue; + std::vector seen(graph.size(), false); + std::vector start_index(parent.size()); + for (int i = 0; i < parent.size(); ++i) { + // Note that because of the way we constructed 'parent', the graph is a + // binary tree. This is not required for the correctness of the algorithm + // here though. + CHECK(graph[i].empty() || graph[i].size() == 2); + if (parent[i] != i) continue; + + // Explore the subtree rooted at node i. + CHECK(!seen[i]); + queue.push_back(i); + while (!queue.empty()) { + const int node = queue.back(); + if (seen[node]) { + queue.pop_back(); + // All the children of node are in the span [start, end) of the + // subset_data vector. + const int start = start_index[node]; + if (new_size - start >= min_subset_size) { + subsets->emplace_back(&(*subset_data)[start], new_size - start); + } + continue; + } + seen[node] = true; + start_index[node] = new_size; + if (node < num_nodes) (*subset_data)[new_size++] = node; + for (const int child : graph[node]) { + if (!seen[child]) queue.push_back(child); + } + } + } + } + + DCHECK_EQ(new_size, num_nodes); +} + // We roughly follow the algorithm described in section 6 of "The Traveling // Salesman Problem, A computational Study", David L. Applegate, Robert E. // Bixby, Vasek Chvatal, William J. Cook. @@ -2423,99 +2526,16 @@ void SeparateSubtourInequalities( arc_by_decreasing_lp_values.end(), std::greater>()); - // We will do a union-find by adding one by one the arc of the lp solution - // in the order above. Every intermediate set during this construction will - // be a candidate for a cut. - // - // In parallel to the union-find, to efficiently reconstruct these sets (at - // most num_nodes), we construct a "decomposition forest" of the different - // connected components. Note that we don't exploit any asymmetric nature of - // the graph here. This is exactly the algo 6.3 in the book above. - int num_components = num_nodes; - std::vector parent(num_nodes); - std::vector root(num_nodes); - for (int i = 0; i < num_nodes; ++i) { - parent[i] = i; - root[i] = i; + std::vector> ordered_arcs; + for (const auto& [score, arc] : arc_by_decreasing_lp_values) { + ordered_arcs.push_back({tails[arc], heads[arc]}); } - auto get_root_and_compress_path = [&root](int node) { - int r = node; - while (root[r] != r) r = root[r]; - while (root[node] != r) { - const int next = root[node]; - root[node] = r; - node = next; - } - return r; - }; - for (const auto& pair : arc_by_decreasing_lp_values) { - if (num_components == 2) break; - const int tail = get_root_and_compress_path(tails[pair.second]); - const int head = get_root_and_compress_path(heads[pair.second]); - if (tail != head) { - // Update the decomposition forest, note that the number of nodes is - // growing. - const int new_node = parent.size(); - parent.push_back(new_node); - parent[head] = new_node; - parent[tail] = new_node; - --num_components; - - // It is important that the union-find representative is the same node. - root.push_back(new_node); - root[head] = new_node; - root[tail] = new_node; - } - } - - // For each node in the decomposition forest, try to add a cut for the set - // formed by the nodes and its children. To do that efficiently, we first - // order the nodes so that for each node in a tree, the set of children forms - // a consecutive span in the pre_order vector. This vector just lists the - // nodes in the "pre-order" graph traversal order. The Spans will point inside - // the pre_order vector, it is why we initialize it once and for all. - int new_size = 0; - std::vector pre_order(num_nodes); + std::vector subset_data; std::vector> subsets; - { - std::vector> graph(parent.size()); - for (int i = 0; i < parent.size(); ++i) { - if (parent[i] != i) graph[parent[i]].push_back(i); - } - std::vector queue; - std::vector seen(graph.size(), false); - std::vector start_index(parent.size()); - for (int i = num_nodes; i < parent.size(); ++i) { - // Note that because of the way we constructed 'parent', the graph is a - // binary tree. This is not required for the correctness of the algorithm - // here though. - CHECK(graph[i].empty() || graph[i].size() == 2); - if (parent[i] != i) continue; - - // Explore the subtree rooted at node i. - CHECK(!seen[i]); - queue.push_back(i); - while (!queue.empty()) { - const int node = queue.back(); - if (seen[node]) { - queue.pop_back(); - // All the children of node are in the span [start, end) of the - // pre_order vector. - const int start = start_index[node]; - if (new_size - start > 1) { - subsets.emplace_back(&pre_order[start], new_size - start); - } - continue; - } - seen[node] = true; - start_index[node] = new_size; - if (node < num_nodes) pre_order[new_size++] = node; - for (const int child : graph[node]) { - if (!seen[child]) queue.push_back(child); - } - } - } - } + GenerateInterestingSubsets(num_nodes, ordered_arcs, + /*min_subset_size=*/2, + /*stop_at_num_components=*/2, &subset_data, + &subsets); // Compute the total demands in order to know the minimum incoming/outgoing // flow. @@ -2525,11 +2545,10 @@ void SeparateSubtourInequalities( } // Process each subsets and add any violated cut. - CHECK_EQ(pre_order.size(), num_nodes); std::vector in_subset(num_nodes, false); for (const absl::Span subset : subsets) { - CHECK_GT(subset.size(), 1); - CHECK_LT(subset.size(), num_nodes); + DCHECK_GT(subset.size(), 1); + DCHECK_LT(subset.size(), num_nodes); // These fields will be left untouched if demands.empty(). bool contain_depot = false; @@ -2661,6 +2680,175 @@ CutGenerator CreateCVRPCutGenerator(int num_nodes, return result; } +// This is really similar to SeparateSubtourInequalities, see the reference +// there. +void SeparateFlowInequalities( + int num_nodes, const std::vector& tails, const std::vector& heads, + const std::vector& arc_capacities, + std::function& in_subset, + IntegerValue* min_incoming_flow, + IntegerValue* min_outgoing_flow)> + get_flows, + const absl::StrongVector& lp_values, + LinearConstraintManager* manager, Model* model) { + // We will collect only the arcs with a positive lp capacity value to speed up + // some computation below. + struct Arc { + int tail; + int head; + double lp_value; + IntegerValue offset; + }; + std::vector relevant_arcs; + + // Often capacities have a coeff > 1. + // We currently exploit this if all coeff have a gcd > 1. + int64_t gcd = 0; + + // Sort the arcs by non-increasing lp_values. + std::vector> arc_by_decreasing_lp_values; + for (int i = 0; i < arc_capacities.size(); ++i) { + const double lp_value = arc_capacities[i].LpValue(lp_values); + if (!arc_capacities[i].IsConstant()) { + gcd = MathUtil::GCD64(gcd, std::abs(arc_capacities[i].coeff.value())); + } + if (lp_value < 1e-6 && arc_capacities[i].constant == 0) continue; + relevant_arcs.push_back( + {tails[i], heads[i], lp_value, arc_capacities[i].constant}); + arc_by_decreasing_lp_values.push_back({lp_value, i}); + } + if (gcd == 0) return; + std::sort(arc_by_decreasing_lp_values.begin(), + arc_by_decreasing_lp_values.end(), + std::greater>()); + + std::vector> ordered_arcs; + for (const auto& [score, arc] : arc_by_decreasing_lp_values) { + if (tails[arc] == -1) continue; + if (heads[arc] == -1) continue; + ordered_arcs.push_back({tails[arc], heads[arc]}); + } + std::vector subset_data; + std::vector> subsets; + GenerateInterestingSubsets(num_nodes, ordered_arcs, + /*min_subset_size=*/1, + /*stop_at_num_components=*/1, &subset_data, + &subsets); + + // Process each subsets and add any violated cut. + std::vector in_subset(num_nodes, false); + for (const absl::Span subset : subsets) { + DCHECK(!subset.empty()); + DCHECK_LE(subset.size(), num_nodes); + + // Initialize "in_subset" and the subset demands. + for (const int n : subset) in_subset[n] = true; + + IntegerValue min_incoming_flow; + IntegerValue min_outgoing_flow; + get_flows(in_subset, &min_incoming_flow, &min_outgoing_flow); + + // We will sum the offset of all incoming/outgoing arc capacities. + // Note that all arcs with a non-zero offset are part of relevant_arcs. + IntegerValue incoming_offset(0); + IntegerValue outgoing_offset(0); + + // Compute the current flow in and out of the subset. + // + // This can take a significant portion of the running time, it is why it is + // faster to do it only on arcs with non-zero lp values which should be in + // linear number rather than the total number of arc which can be quadratic. + double lp_outgoing_flow = 0.0; + double lp_incoming_flow = 0.0; + for (const auto arc : relevant_arcs) { + if (arc.tail != -1 && in_subset[arc.tail]) { + if (arc.head == -1 || !in_subset[arc.head]) { + incoming_offset += arc.offset; + lp_outgoing_flow += arc.lp_value; + } + } else { + if (arc.head != -1 && in_subset[arc.head]) { + outgoing_offset += arc.offset; + lp_incoming_flow += arc.lp_value; + } + } + } + + // If the gcd is greater than one, because all variables are integer we + // can round the flow lower bound to the next multiple of the gcd. + // + // TODO(user): Alternatively, try MIR heuristics if the coefficients in + // the capacities are not all the same. + if (gcd > 1) { + const IntegerValue test_incoming = min_incoming_flow - incoming_offset; + const IntegerValue new_incoming = + CeilRatio(test_incoming, IntegerValue(gcd)) * IntegerValue(gcd); + const IntegerValue incoming_delta = new_incoming - test_incoming; + if (incoming_delta > 0) min_incoming_flow += incoming_delta; + } + if (gcd > 1) { + const IntegerValue test_outgoing = min_outgoing_flow - outgoing_offset; + const IntegerValue new_outgoing = + CeilRatio(test_outgoing, IntegerValue(gcd)) * IntegerValue(gcd); + const IntegerValue outgoing_delta = new_outgoing - test_outgoing; + if (outgoing_delta > 0) min_outgoing_flow += outgoing_delta; + } + + if (lp_incoming_flow < ToDouble(min_incoming_flow) - 1e-6) { + VLOG(2) << "INCOMING CUT " << lp_incoming_flow + << " >= " << min_incoming_flow << " size " << subset.size() + << " offset " << incoming_offset << " gcd " << gcd; + LinearConstraintBuilder cut(model, min_incoming_flow, kMaxIntegerValue); + for (int i = 0; i < tails.size(); ++i) { + if ((tails[i] == -1 || !in_subset[tails[i]]) && + (heads[i] != -1 && in_subset[heads[i]])) { + cut.AddTerm(arc_capacities[i], 1.0); + } + } + manager->AddCut(cut.Build(), "IncomingFlow", lp_values); + } + + if (lp_outgoing_flow < ToDouble(min_outgoing_flow) - 1e-6) { + VLOG(2) << "OUGOING CUT " << lp_outgoing_flow + << " >= " << min_outgoing_flow << " size " << subset.size() + << " offset " << outgoing_offset << " gcd " << gcd; + LinearConstraintBuilder cut(model, min_outgoing_flow, kMaxIntegerValue); + for (int i = 0; i < tails.size(); ++i) { + if ((tails[i] != -1 && in_subset[tails[i]]) && + (heads[i] == -1 || !in_subset[heads[i]])) { + cut.AddTerm(arc_capacities[i], 1.0); + } + } + manager->AddCut(cut.Build(), "OutgoingFlow", lp_values); + } + + // Sparse clean up. + for (const int n : subset) in_subset[n] = false; + } +} + +CutGenerator CreateFlowCutGenerator( + int num_nodes, const std::vector& tails, const std::vector& heads, + const std::vector& arc_capacities, + std::function& in_subset, + IntegerValue* min_incoming_flow, + IntegerValue* min_outgoing_flow)> + get_flows, + Model* model) { + CutGenerator result; + for (const AffineExpression expr : arc_capacities) { + if (!expr.IsConstant()) result.vars.push_back(expr.var); + } + result.generate_cuts = + [=](const absl::StrongVector& lp_values, + LinearConstraintManager* manager) { + SeparateFlowInequalities(num_nodes, tails, heads, arc_capacities, + get_flows, lp_values, manager, model); + return true; + }; + return result; +} + std::function LinearProgrammingConstraint::HeuristicLpMostInfeasibleBinary(Model* model) { // Gather all 0-1 variables that appear in this LP. diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 97b2aba65a..3b7a86664a 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -569,6 +569,35 @@ class LinearProgrammingDispatcher class LinearProgrammingConstraintCollection : public std::vector {}; +// TODO(user): Move the cut "graph" based cut generator out of this class, there +// is no reason to keep them here. + +// Given a graph with nodes in [0, num_nodes) and a set of arcs (the order is +// important), this will: +// - Start with each nodes in separate "subsets". +// - Consider the arc in order, and each time one connects two separate +// subsets, merge the two subsets into a new one. +// - Stops when there is only 2 subset left. +// - Output all subsets generated this way (at most 2 * num_nodes). The +// subsets spans will point in the subset_data vector (which will be of size +// exactly num_nodes). +// +// Only subsets of size >= min_subset_size will be returned. This is mainly here +// to exclude subsets of size 1. +// +// This is an heuristic to generate interesting cuts for TSP or other graph +// based constraints. We roughly follow the algorithm described in section 6 of +// "The Traveling Salesman Problem, A computational Study", David L. Applegate, +// Robert E. Bixby, Vasek Chvatal, William J. Cook. +// +// Note that this is mainly a "symmetric" case algo, but it does still work for +// the asymmetric case. +void GenerateInterestingSubsets(int num_nodes, + const std::vector>& arcs, + int min_subset_size, int stop_at_num_components, + std::vector* subset_data, + std::vector>* subsets); + // Cut generator for the circuit constraint, where in any feasible solution, the // arcs that are present (variable at 1) must form a circuit through all the // nodes of the graph. Self arc are forbidden in this case. @@ -590,6 +619,26 @@ CutGenerator CreateCVRPCutGenerator(int num_nodes, const std::vector& literals, const std::vector& demands, int64_t capacity, Model* model); + +// Try to find a subset where the current LP capacity of the outgoing or +// incoming arc is not enough to satisfy the demands. +// +// We support the special value -1 for tail or head that means that the arc +// comes from (or is going to) outside the nodes in [0, num_nodes). Such arc +// must still have a capacity assigned to it. +// +// TODO(user): Support general linear expression for capacities. +// TODO(user): Some model applies the same capacity to both an arc and its +// reverse. Also support this case. +CutGenerator CreateFlowCutGenerator( + int num_nodes, const std::vector& tails, const std::vector& heads, + const std::vector& arc_capacities, + std::function& in_subset, + IntegerValue* min_incoming_flow, + IntegerValue* min_outgoing_flow)> + get_flows, + Model* model); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/scheduling_constraints.cc b/ortools/sat/scheduling_constraints.cc index ab7aa38473..c6a046746e 100644 --- a/ortools/sat/scheduling_constraints.cc +++ b/ortools/sat/scheduling_constraints.cc @@ -343,6 +343,7 @@ std::function SpanOfIntervals( std::vector starts; std::vector ends; std::vector clause; + absl::flat_hash_set added_literals; bool at_least_one_interval_is_present = false; const Literal true_literal = model->GetOrCreate()->GetTrueLiteral(); @@ -353,7 +354,9 @@ std::function SpanOfIntervals( if (repository->IsOptional(interval)) { const Literal task_lit = repository->PresenceLiteral(interval); presence_literals.push_back(task_lit); - clause.push_back(task_lit); + if (added_literals.insert(task_lit.Index()).second) { + clause.push_back(task_lit); + } if (repository->IsOptional(span)) { // task is present => target is present. @@ -372,7 +375,10 @@ std::function SpanOfIntervals( if (!at_least_one_interval_is_present) { // enforcement_literal is true => one of the task is present. if (repository->IsOptional(span)) { - clause.push_back(repository->PresenceLiteral(span).Negated()); + const Literal span_literal = repository->PresenceLiteral(span); + if (added_literals.insert(span_literal.Negated().Index()).second) { + clause.push_back(span_literal.Negated()); + } } sat_solver->AddProblemClause(clause); }