diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index 1fca363bfb..f198049d77 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -722,7 +722,10 @@ cc_library( "//ortools/util:sorted_interval_list", "//ortools/util:strong_integers", "//ortools/util:time_limit", + "@com_google_absl//absl/container:btree", + "@com_google_absl//absl/container:flat_hash_set", "@com_google_absl//absl/container:inlined_vector", + "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", "@com_google_absl//absl/strings", "@com_google_absl//absl/types:span", @@ -792,11 +795,9 @@ cc_library( ":sat_parameters_cc_proto", ":util", "//ortools/base", - "//ortools/base:hash", "//ortools/base:stl_util", "//ortools/base:strong_vector", "//ortools/base:timer", - "//ortools/base:types", "//ortools/graph:strongly_connected_components", "//ortools/util:bitset", "//ortools/util:stats", @@ -932,10 +933,7 @@ cc_library( ":sat_solver", "//ortools/base", "//ortools/base:cleanup", - "//ortools/base:hash", "//ortools/base:strong_vector", - "//ortools/base:types", - "//ortools/graph:iterators", "//ortools/util:bitset", "//ortools/util:rev", "//ortools/util:saturated_arithmetic", @@ -958,6 +956,7 @@ cc_library( srcs = ["integer_search.cc"], hdrs = ["integer_search.h"], deps = [ + ":clause", ":cp_model_cc_proto", ":cp_model_mapping", ":implied_bounds", @@ -984,7 +983,7 @@ cc_library( "@com_google_absl//absl/meta:type_traits", "@com_google_absl//absl/random:distributions", "@com_google_absl//absl/strings", - "@com_google_absl//absl/time", + "@com_google_absl//absl/types:span", ], ) @@ -1501,6 +1500,7 @@ cc_library( "//ortools/util:strong_integers", "@com_google_absl//absl/container:inlined_vector", "@com_google_absl//absl/log:check", + "@com_google_absl//absl/strings", "@com_google_absl//absl/types:span", ], ) @@ -1789,12 +1789,13 @@ cc_library( ":intervals", "//ortools/base", "//ortools/base:stl_util", - "//ortools/graph:connected_components", "//ortools/util:integer_pq", "//ortools/util:strong_integers", "@com_google_absl//absl/container:flat_hash_set", + "@com_google_absl//absl/container:inlined_vector", "@com_google_absl//absl/log:check", "@com_google_absl//absl/random:bit_gen_ref", + "@com_google_absl//absl/random:distributions", "@com_google_absl//absl/strings:str_format", "@com_google_absl//absl/types:span", ], @@ -1813,13 +1814,15 @@ cc_library( ":intervals", ":linear_constraint", ":model", - ":precedences", ":sat_base", ":sat_parameters_cc_proto", + ":synchronization", ":timetable", + ":util", "//ortools/util:saturated_arithmetic", "//ortools/util:strong_integers", "@com_google_absl//absl/container:flat_hash_set", + "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", "@com_google_absl//absl/types:span", ], @@ -2200,15 +2203,3 @@ cc_library( "@com_google_absl//absl/types:span", ], ) - -cc_library( - name = "cp_model_objective", - srcs = ["cp_model_objective.cc"], - hdrs = ["cp_model_objective.h"], - deps = [ - ":cp_model_cc_proto", - ":cp_model_utils", - "//ortools/util:sorted_interval_list", - "@com_google_absl//absl/log:check", - ], -) diff --git a/ortools/sat/circuit.h b/ortools/sat/circuit.h index 5d57907ad3..c7a98c4c57 100644 --- a/ortools/sat/circuit.h +++ b/ortools/sat/circuit.h @@ -244,7 +244,7 @@ int ReindexArcs(IntContainer* tails, IntContainer* heads, // ============================================================================ // This just wraps CircuitPropagator. See the comment there to see what this -// does. Note that any nodes with no outoing or no incoming arc will cause the +// does. Note that any nodes with no outgoing or no incoming arc will cause the // problem to be UNSAT. One can call ReindexArcs() first to ignore such nodes. std::function SubcircuitConstraint( int num_nodes, const std::vector& tails, const std::vector& heads, diff --git a/ortools/sat/clause.cc b/ortools/sat/clause.cc index c6cc9ae7d5..9e74d8d27e 100644 --- a/ortools/sat/clause.cc +++ b/ortools/sat/clause.cc @@ -39,11 +39,11 @@ #include "ortools/sat/inclusion.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" +#include "ortools/sat/util.h" #include "ortools/util/bitset.h" #include "ortools/util/stats.h" #include "ortools/util/strong_integers.h" #include "ortools/util/time_limit.h" - namespace operations_research { namespace sat { @@ -207,7 +207,7 @@ bool LiteralWatchers::Propagate(Trail* trail) { return true; } -absl::Span LiteralWatchers::Reason(const Trail& trail, +absl::Span LiteralWatchers::Reason(const Trail& /*trail*/, int trail_index) const { return reasons_[trail_index]->PropagationReason(); } @@ -291,8 +291,8 @@ bool LiteralWatchers::AttachAndPropagate(SatClause* clause, Trail* trail) { void LiteralWatchers::Attach(SatClause* clause, Trail* trail) { Literal* literals = clause->literals(); - CHECK(!trail->Assignment().LiteralIsAssigned(literals[0])); - CHECK(!trail->Assignment().LiteralIsAssigned(literals[1])); + DCHECK(!trail->Assignment().LiteralIsAssigned(literals[0])); + DCHECK(!trail->Assignment().LiteralIsAssigned(literals[1])); ++num_watched_clauses_; AttachOnFalse(literals[0], literals[1], clause); @@ -321,7 +321,7 @@ void LiteralWatchers::Detach(SatClause* clause) { for (const Literal l : {clause->FirstLiteral(), clause->SecondLiteral()}) { needs_cleaning_.Clear(l); RemoveIf(&(watchers_on_false_[l]), [](const Watcher& watcher) { - return !watcher.clause->IsAttached(); + return watcher.clause->IsRemoved(); }); } } @@ -347,7 +347,7 @@ void LiteralWatchers::AttachAllClauses() { DeleteRemovedClauses(); for (SatClause* clause : clauses_) { ++num_watched_clauses_; - CHECK_GE(clause->size(), 2); + DCHECK_GE(clause->size(), 2); AttachOnFalse(clause->FirstLiteral(), clause->SecondLiteral(), clause); AttachOnFalse(clause->SecondLiteral(), clause->FirstLiteral(), clause); } @@ -355,7 +355,7 @@ void LiteralWatchers::AttachAllClauses() { // This one do not need the clause to be detached. bool LiteralWatchers::InprocessingFixLiteral(Literal true_literal) { - CHECK_EQ(trail_->CurrentDecisionLevel(), 0); + DCHECK_EQ(trail_->CurrentDecisionLevel(), 0); if (drat_proof_handler_ != nullptr) { drat_proof_handler_->AddClause({true_literal}); } @@ -374,7 +374,7 @@ bool LiteralWatchers::InprocessingFixLiteral(Literal true_literal) { // TODO(user): We could do something slower if the clauses are attached like // we do for InprocessingRewriteClause(). void LiteralWatchers::InprocessingRemoveClause(SatClause* clause) { - CHECK(!all_clauses_are_attached_); + DCHECK(!all_clauses_are_attached_); if (drat_proof_handler_ != nullptr) { drat_proof_handler_->DeleteClause(clause->AsSpan()); } @@ -388,7 +388,7 @@ bool LiteralWatchers::InprocessingRewriteClause( if (DEBUG_MODE) { for (const Literal l : new_clause) { - CHECK(!trail_->Assignment().LiteralIsAssigned(l)); + DCHECK(!trail_->Assignment().LiteralIsAssigned(l)); } } @@ -399,7 +399,7 @@ bool LiteralWatchers::InprocessingRewriteClause( } if (new_clause.size() == 2) { - implication_graph_->AddBinaryClause(new_clause[0], new_clause[1]); + CHECK(implication_graph_->AddBinaryClause(new_clause[0], new_clause[1])); InprocessingRemoveClause(clause); return true; } @@ -418,25 +418,25 @@ bool LiteralWatchers::InprocessingRewriteClause( for (const Literal l : {clause->FirstLiteral(), clause->SecondLiteral()}) { needs_cleaning_.Clear(l); RemoveIf(&(watchers_on_false_[l]), [](const Watcher& watcher) { - return !watcher.clause->IsAttached(); + return watcher.clause->IsRemoved(); }); } } clause->Rewrite(new_clause); - // And we re-attach it. + // And we reattach it. if (all_clauses_are_attached_) Attach(clause, trail_); return true; } SatClause* LiteralWatchers::InprocessingAddClause( absl::Span new_clause) { - CHECK(!new_clause.empty()); - CHECK(!all_clauses_are_attached_); + DCHECK(!new_clause.empty()); + DCHECK(!all_clauses_are_attached_); if (DEBUG_MODE) { for (const Literal l : new_clause) { - CHECK(!trail_->Assignment().LiteralIsAssigned(l)); + DCHECK(!trail_->Assignment().LiteralIsAssigned(l)); } } @@ -461,7 +461,7 @@ void LiteralWatchers::CleanUpWatchers() { for (const LiteralIndex index : needs_cleaning_.PositionsSetAtLeastOnce()) { DCHECK(needs_cleaning_[index]); RemoveIf(&(watchers_on_false_[index]), [](const Watcher& watcher) { - return !watcher.clause->IsAttached(); + return watcher.clause->IsRemoved(); }); needs_cleaning_.Clear(index); } @@ -479,13 +479,13 @@ void LiteralWatchers::DeleteRemovedClauses() { to_minimize_index_ = std::stable_partition(clauses_.begin(), clauses_.begin() + to_minimize_index_, - [](SatClause* a) { return a->IsAttached(); }) - + [](SatClause* a) { return !a->IsRemoved(); }) - clauses_.begin(); // Do the proper deletion. std::vector::iterator iter = std::stable_partition(clauses_.begin(), clauses_.end(), - [](SatClause* a) { return a->IsAttached(); }); + [](SatClause* a) { return !a->IsRemoved(); }); gtl::STLDeleteContainerPointers(iter, clauses_.end()); clauses_.erase(iter, clauses_.end()); } @@ -521,6 +521,11 @@ bool BinaryImplicationGraph::AddBinaryClause(Literal a, Literal b) { drat_proof_handler_->AddClause({a, b}); } + if (is_redundant_[a.Index()]) a = Literal(representative_of_[a.Index()]); + if (is_redundant_[b.Index()]) b = Literal(representative_of_[b.Index()]); + if (a == b.Negated()) return true; + if (a == b) return FixLiteral(a); + DCHECK(!is_removed_[a]); DCHECK(!is_removed_[b]); estimated_sizes_[a.NegatedIndex()]++; @@ -532,8 +537,8 @@ bool BinaryImplicationGraph::AddBinaryClause(Literal a, Literal b) { const auto& assignment = trail_->Assignment(); if (trail_->CurrentDecisionLevel() == 0) { - CHECK(!assignment.LiteralIsAssigned(a)); - CHECK(!assignment.LiteralIsAssigned(b)); + DCHECK(!assignment.LiteralIsAssigned(a)); + DCHECK(!assignment.LiteralIsAssigned(b)); } else { if (assignment.LiteralIsFalse(a)) { if (assignment.LiteralIsAssigned(b)) { @@ -553,9 +558,9 @@ bool BinaryImplicationGraph::AddBinaryClause(Literal a, Literal b) { return true; } -bool BinaryImplicationGraph::AddAtMostOne(absl::Span at_most_one, - int expansion_size) { - CHECK_EQ(trail_->CurrentDecisionLevel(), 0); +bool BinaryImplicationGraph::AddAtMostOne( + absl::Span at_most_one) { + DCHECK_EQ(trail_->CurrentDecisionLevel(), 0); if (at_most_one.size() <= 1) return true; // Temporarily copy the at_most_one constraint at the end of @@ -567,12 +572,15 @@ bool BinaryImplicationGraph::AddAtMostOne(absl::Span at_most_one, at_most_one_buffer_.push_back(Literal(kNoLiteralIndex)); is_dag_ = false; - return CleanUpAndAddAtMostOnes(base_index, expansion_size); + return CleanUpAndAddAtMostOnes(base_index); } -// TODO(user): remove duplication with -// LiteralWatchers::InprocessingFixLiteral(); +// TODO(user): remove dupl with LiteralWatchers::InprocessingFixLiteral(). +// +// Note that we currently do not support calling this at a positive level since +// we might loose the fixing on backtrack otherwise. bool BinaryImplicationGraph::FixLiteral(Literal true_literal) { + CHECK_EQ(trail_->CurrentDecisionLevel(), 0); if (trail_->Assignment().LiteralIsTrue(true_literal)) return true; if (trail_->Assignment().LiteralIsFalse(true_literal)) return false; @@ -587,8 +595,7 @@ bool BinaryImplicationGraph::FixLiteral(Literal true_literal) { // This works by doing a linear scan on the at_most_one_buffer_ and // cleaning/copying the at most ones on the fly to the beginning of the same // buffer. -bool BinaryImplicationGraph::CleanUpAndAddAtMostOnes(const int base_index, - int expansion_size) { +bool BinaryImplicationGraph::CleanUpAndAddAtMostOnes(int base_index) { const VariablesAssignment& assignment = trail_->Assignment(); int local_end = base_index; const int buffer_size = at_most_one_buffer_.size(); @@ -598,6 +605,17 @@ bool BinaryImplicationGraph::CleanUpAndAddAtMostOnes(const int base_index, // Process a new at most one. // It will be copied into buffer[local_start, local_end]. const int local_start = local_end; + + // Update the iterator now. Even if the current at_most_one is reduced away, + // local_start will still point to the next one, or to the end of the + // buffer. + if (i == at_most_one_iterator_) { + at_most_one_iterator_ = local_start; + DCHECK(at_most_one_iterator_ == 0 || + at_most_one_buffer_[at_most_one_iterator_ - 1].Index() == + kNoLiteralIndex); + } + bool set_all_left_to_false = false; for (;; ++i) { const Literal l = at_most_one_buffer_[i]; @@ -679,8 +697,7 @@ bool BinaryImplicationGraph::CleanUpAndAddAtMostOnes(const int base_index, // We expand small sizes into implications. // TODO(user): Investigate what the best threshold is. - if (at_most_one.size() < expansion_size) { - // Note that his automatically skip size 0 and 1. + if (at_most_one.size() <= std::max(2, at_most_one_max_expansion_size_)) { for (const Literal a : at_most_one) { for (const Literal b : at_most_one) { if (a == b) continue; @@ -699,7 +716,7 @@ bool BinaryImplicationGraph::CleanUpAndAddAtMostOnes(const int base_index, if (l.Index() >= at_most_ones_.size()) { at_most_ones_.resize(l.Index().value() + 1); } - CHECK(!is_redundant_[l]); + DCHECK(!is_redundant_[l]); at_most_ones_[l].push_back(local_start); } @@ -756,7 +773,7 @@ bool BinaryImplicationGraph::PropagateOnTrue(Literal true_literal, ++num_inspections_; if (literal == true_literal) { if (DEBUG_MODE) { - CHECK(!seen); + DCHECK(!seen); seen = true; } continue; @@ -795,7 +812,7 @@ bool BinaryImplicationGraph::Propagate(Trail* trail) { } absl::Span BinaryImplicationGraph::Reason( - const Trail& trail, int trail_index) const { + const Trail& /*trail*/, int trail_index) const { return {&reasons_[trail_index], 1}; } @@ -887,7 +904,7 @@ void BinaryImplicationGraph::MinimizeConflictFirst( const Trail& trail, std::vector* conflict, SparseBitset* marked) { SCOPED_TIME_STAT(&stats_); - CHECK(!conflict->empty()); + DCHECK(!conflict->empty()); is_marked_.ClearAndResize(LiteralIndex(implications_.size())); MarkDescendants(conflict->front().Negated()); for (const LiteralIndex i : is_marked_.PositionsSetAtLeastOnce()) { @@ -902,7 +919,7 @@ void BinaryImplicationGraph::MinimizeConflictFirst( // computation to remove redundant implication in the implication list of the // first UIP conflict. void BinaryImplicationGraph::MinimizeConflictFirstWithTransitiveReduction( - const Trail& trail, std::vector* conflict, + const Trail& /*trail*/, std::vector* conflict, absl::BitGenRef random) { SCOPED_TIME_STAT(&stats_); const LiteralIndex root_literal_index = conflict->front().NegatedIndex(); @@ -1016,7 +1033,7 @@ void BinaryImplicationGraph::MinimizeConflictExperimental( void BinaryImplicationGraph::RemoveFixedVariables() { SCOPED_TIME_STAT(&stats_); - CHECK_EQ(trail_->CurrentDecisionLevel(), 0); + DCHECK_EQ(trail_->CurrentDecisionLevel(), 0); if (IsEmpty()) return; // Nothing to do if nothing changed since last call. @@ -1033,7 +1050,7 @@ void BinaryImplicationGraph::RemoveFixedVariables() { // The code assumes that everything is already propagated. // Otherwise we will remove implications that didn't propagate yet! for (const Literal lit : implications_[true_literal]) { - CHECK(trail_->Assignment().LiteralIsTrue(lit)); + DCHECK(trail_->Assignment().LiteralIsTrue(lit)); } } @@ -1074,6 +1091,7 @@ void BinaryImplicationGraph::RemoveFixedVariables() { // this shouldn't change the correctness of the code. at_most_ones_.clear(); CleanUpAndAddAtMostOnes(/*base_index=*/0); + DCHECK(InvariantsAreOk()); } class SccGraph { @@ -1119,7 +1137,7 @@ class SccGraph { if (at_most_one_already_explored_[start]) { // We never expand a node twice. const int first_node = previous_node_to_explore_at_most_one_[start]; - CHECK_NE(node, first_node); + DCHECK_NE(node, first_node); if (finder_.NodeIsInCurrentDfsPath(first_node)) { // If the first node is not settled, then we do explore the @@ -1374,7 +1392,7 @@ bool BinaryImplicationGraph::DetectEquivalences(bool log_info) { // TODO(user): Track which literal have new implications, and only process // the antecedants of these. bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { - CHECK_EQ(trail_->CurrentDecisionLevel(), 0); + DCHECK_EQ(trail_->CurrentDecisionLevel(), 0); if (!DetectEquivalences()) return false; // TODO(user): the situation with fixed variable is not really "clean". @@ -1463,7 +1481,7 @@ bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { // We have a DAG, so direct_child could only be marked first. is_marked_.Clear(direct_child); } - CHECK(!is_marked_[root]) + DCHECK(!is_marked_[root]) << "DetectEquivalences() should have removed cycles!"; is_marked_.Set(root); @@ -1490,10 +1508,10 @@ bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { const LiteralIndex i = marked_positions[marked_index]; if (is_marked_[Literal(i).NegatedIndex()]) { // We tested that at the beginning. - CHECK(!trail_->Assignment().LiteralIsAssigned(Literal(root))); + DCHECK(!trail_->Assignment().LiteralIsAssigned(Literal(root))); // We propagate right away the binary implications so that we do not - // need to consider all antecedants of root in the transitive + // need to consider all antecedents of root in the transitive // reduction. ++num_fixed; if (!FixLiteral(Literal(root).Negated())) return false; @@ -1515,7 +1533,7 @@ bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { direct_implications[new_size++] = l; } else { tmp_removed_.push_back({Literal(root), l}); - CHECK(!is_redundant_[l]); + DCHECK(!is_redundant_[l]); } } const int diff = direct_implications.size() - new_size; @@ -1534,11 +1552,8 @@ bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { is_marked_.ClearAndResize(size); // If we aborted early, we might no longer have both a=>b and not(b)=>not(a). - // This is not desirable has some algo relies on this invariant. That code fix - // this. - if (num_fixed > 0) { - RemoveFixedVariables(); - } + // This is not desirable has some algo relies on this invariant. We fix this + // here. if (aborted) { absl::flat_hash_set> removed; for (const auto [a, b] : tmp_removed_) { @@ -1555,6 +1570,9 @@ bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { implication.resize(new_size); } } + if (num_fixed > 0) { + RemoveFixedVariables(); + } DCHECK(InvariantsAreOk()); gtl::STLClearObject(&tmp_removed_); @@ -1666,7 +1684,7 @@ bool BinaryImplicationGraph::TransformIntoMaxCliques( // Special case for clique of size 2, we don't expand them if they // are included in an already added clique. if (clique.size() == 2) { - CHECK_NE(clique[0], clique[1]); + DCHECK_NE(clique[0], clique[1]); const int dense_index = ElementInIntersectionOrMinusOne( max_cliques_containing[clique[0]], max_cliques_containing[clique[1]]); if (dense_index >= 0) { @@ -2082,7 +2100,7 @@ std::vector BinaryImplicationGraph::ExpandAtMostOne( // TODO(user): Mark fixed variable as is_removed_ for faster iteration? const std::vector& BinaryImplicationGraph::DirectImplications( Literal literal) { - CHECK(!is_removed_[literal]); + DCHECK(!is_removed_[literal]); // Clear old state. for (const Literal l : direct_implications_) { @@ -2092,7 +2110,7 @@ const std::vector& BinaryImplicationGraph::DirectImplications( // Fill new state. const VariablesAssignment& assignment = trail_->Assignment(); - CHECK(!assignment.LiteralIsAssigned(literal)); + DCHECK(!assignment.LiteralIsAssigned(literal)); for (const Literal l : implications_[literal]) { if (l == literal) continue; if (assignment.LiteralIsAssigned(l)) continue; @@ -2103,7 +2121,7 @@ const std::vector& BinaryImplicationGraph::DirectImplications( } if (literal.Index() < at_most_ones_.size()) { if (is_redundant_[literal]) { - CHECK(at_most_ones_[literal].empty()); + DCHECK(at_most_ones_[literal].empty()); } for (const int start : at_most_ones_[literal]) { for (int i = start;; ++i) { @@ -2125,7 +2143,7 @@ const std::vector& BinaryImplicationGraph::DirectImplications( bool BinaryImplicationGraph::FindFailedLiteralAroundVar(BooleanVariable var, bool* is_unsat) { const int saved_index = propagation_trail_index_; - CHECK_EQ(propagation_trail_index_, trail_->Index()); // Propagation done. + DCHECK_EQ(propagation_trail_index_, trail_->Index()); // Propagation done. const VariablesAssignment& assignment = trail_->Assignment(); if (assignment.VariableIsAssigned(var)) return false; @@ -2158,7 +2176,7 @@ int64_t BinaryImplicationGraph::NumImplicationOnVariableRemoval( result += s1; // We should have dealt with that in FindFailedLiteralAroundVar(). - CHECK(!in_direct_implications_[l]); + DCHECK(!in_direct_implications_[l]); // l => literal => l: equivalent variable! if (in_direct_implications_[l.NegatedIndex()]) result--; @@ -2170,10 +2188,14 @@ int64_t BinaryImplicationGraph::NumImplicationOnVariableRemoval( void BinaryImplicationGraph::RemoveBooleanVariable( BooleanVariable var, std::deque>* postsolve_clauses) { const Literal literal(var, true); + DCHECK(!is_removed_[literal.Index()]); + DCHECK(!is_redundant_[literal.Index()]); + direct_implications_of_negated_literal_ = DirectImplications(literal.Negated()); for (const Literal b : DirectImplications(literal)) { if (is_removed_[b]) continue; + DCHECK(!is_redundant_[b]); estimated_sizes_[b.NegatedIndex()]--; for (const Literal a_negated : direct_implications_of_negated_literal_) { if (a_negated.Negated() == b) continue; @@ -2183,6 +2205,7 @@ void BinaryImplicationGraph::RemoveBooleanVariable( } for (const Literal a_negated : direct_implications_of_negated_literal_) { if (is_removed_[a_negated]) continue; + DCHECK(!is_redundant_[a_negated]); estimated_sizes_[a_negated.NegatedIndex()]--; } @@ -2213,17 +2236,50 @@ void BinaryImplicationGraph::RemoveBooleanVariable( } } -void BinaryImplicationGraph::CleanupAllRemovedVariables() { +void BinaryImplicationGraph::RemoveAllRedundantVariables( + std::deque>* postsolve_clauses) { for (LiteralIndex a(0); a < implications_.size(); ++a) { - if (is_removed_[a]) { - DCHECK(implications_[a].empty()); + if (is_redundant_[a] && !is_removed_[a]) { + postsolve_clauses->push_back( + {Literal(a), Literal(RepresentativeOf(Literal(a))).Negated()}); + is_removed_[a] = true; + gtl::STLClearObject(&(implications_[a])); continue; } int new_size = 0; auto& implication = implications_[a]; for (const Literal l : implication) { - if (!is_removed_[l]) implication[new_size++] = l; + if (!is_redundant_[l]) { + implication[new_size++] = l; + } + } + implication.resize(new_size); + } +} + +void BinaryImplicationGraph::CleanupAllRemovedAndFixedVariables() { + const VariablesAssignment& assignment = trail_->Assignment(); + for (LiteralIndex a(0); a < implications_.size(); ++a) { + if (is_removed_[a] || assignment.LiteralIsAssigned(Literal(a))) { + if (DEBUG_MODE && assignment.LiteralIsTrue(Literal(a))) { + // The code assumes that everything is already propagated. + // Otherwise we will remove implications that didn't propagate yet! + for (const Literal lit : implications_[a]) { + DCHECK(trail_->Assignment().LiteralIsTrue(lit)); + } + } + + gtl::STLClearObject(&(implications_[a])); + continue; + } + + int new_size = 0; + auto& implication = implications_[a]; + for (const Literal l : implication) { + if (!is_removed_[l] && !assignment.LiteralIsTrue(l)) { + implication[new_size++] = l; + } } implication.resize(new_size); } @@ -2231,6 +2287,8 @@ void BinaryImplicationGraph::CleanupAllRemovedVariables() { // Clean-up at most ones. at_most_ones_.clear(); CleanUpAndAddAtMostOnes(/*base_index=*/0); + + // Note that to please the invariant() we also removed fixed literal above. DCHECK(InvariantsAreOk()); } @@ -2248,13 +2306,16 @@ bool BinaryImplicationGraph::InvariantsAreOk() { } continue; } + if (is_removed_[a_index]) { + if (!implications_[a_index].empty()) { + LOG(ERROR) << "Removed literal has non-cleared implications"; + return false; + } + continue; + } if (is_redundant_[a_index]) { ++num_redundant; - - // TODO(user): nothing really prevent other parts of the solver to re-add - // implication containing redundant literal. It should be up to us to - // filter them and replace by the representative... - if (false && implications_[a_index].size() != 1) { + if (implications_[a_index].size() != 1) { LOG(ERROR) << "Redundant literal should only point to its representative " << Literal(a_index) << " => " << implications_[a_index]; @@ -2280,20 +2341,30 @@ bool BinaryImplicationGraph::InvariantsAreOk() { for (LiteralIndex a_index(0); a_index < implications_.size(); ++a_index) { const LiteralIndex not_a_index = Literal(a_index).NegatedIndex(); for (const Literal b : implications_[a_index]) { + if (is_removed_[b]) { + LOG(ERROR) << "A removed literal still appear! " << Literal(a_index) + << " => " << b; + return false; + } + if (!seen.contains({b.NegatedIndex(), not_a_index})) { LOG(ERROR) << "We have " << Literal(a_index) << " => " << b << " but not the reverse implication!"; - LOG(ERROR) << is_redundant_[a_index] << " " - << trail_->Assignment().LiteralIsAssigned(Literal(a_index)); - LOG(ERROR) << is_redundant_[b] << " " - << trail_->Assignment().LiteralIsAssigned(b); + LOG(ERROR) << "redundant[a]: " << is_redundant_[a_index] + << " assigned[a]: " + << trail_->Assignment().LiteralIsAssigned(Literal(a_index)) + << " removed[a]: " << is_removed_[a_index] + << " redundant[b]: " << is_redundant_[b] << " assigned[b]: " + << trail_->Assignment().LiteralIsAssigned(b) + << " removed[b]: " << is_removed_[b]; + return false; } // Test that if we have a dag, our topo order is correct. if (is_dag_ && !is_redundant_[b] && !is_redundant_[a_index]) { - CHECK_NE(lit_to_order[b], -1); - CHECK_LE(lit_to_order[b], lit_to_order[a_index]); + DCHECK_NE(lit_to_order[b], -1); + DCHECK_LE(lit_to_order[b], lit_to_order[a_index]); } } } @@ -2313,6 +2384,11 @@ bool BinaryImplicationGraph::InvariantsAreOk() { start = i + 1; continue; } + if (is_removed_[at_most_one_buffer_[i]]) { + LOG(ERROR) << "A removed literal still appear in an amo " + << at_most_one_buffer_[i]; + return false; + } if (!lit_to_start.contains({at_most_one_buffer_[i], start})) { return false; } @@ -2321,11 +2397,30 @@ bool BinaryImplicationGraph::InvariantsAreOk() { return true; } +absl::Span BinaryImplicationGraph::NextAtMostOne() { + if (at_most_one_iterator_ >= at_most_one_buffer_.size()) { + return absl::Span(); + } + const int local_start = at_most_one_iterator_; + DCHECK_NE(at_most_one_buffer_[local_start].Index(), kNoLiteralIndex); + DCHECK(at_most_one_iterator_ == 0 || + at_most_one_iterator_ >= at_most_one_buffer_.size() || + at_most_one_buffer_[at_most_one_iterator_ - 1].Index() == + kNoLiteralIndex); + int local_end = at_most_one_iterator_ + 1; + while (at_most_one_buffer_[local_end].Index() != kNoLiteralIndex) { + ++local_end; + } + at_most_one_iterator_ = local_end + 1; + return absl::MakeSpan(at_most_one_buffer_.data() + local_start, + local_end - local_start); +} + // ----- SatClause ----- // static SatClause* SatClause::Create(absl::Span literals) { - CHECK_GE(literals.size(), 2); + DCHECK_GE(literals.size(), 2); SatClause* clause = reinterpret_cast( ::operator new(sizeof(SatClause) + literals.size() * sizeof(Literal))); clause->size_ = literals.size(); @@ -2339,7 +2434,7 @@ SatClause* SatClause::Create(absl::Span literals) { // any of the watched literal is assigned, then the clause is necessarily true. bool SatClause::RemoveFixedLiteralsAndTestIfTrue( const VariablesAssignment& assignment) { - DCHECK(IsAttached()); + DCHECK(!IsRemoved()); if (assignment.VariableIsAssigned(literals_[0].Variable()) || assignment.VariableIsAssigned(literals_[1].Variable())) { DCHECK(IsSatisfied(assignment)); diff --git a/ortools/sat/clause.h b/ortools/sat/clause.h index 2d13215553..2a12a1c6e8 100644 --- a/ortools/sat/clause.h +++ b/ortools/sat/clause.h @@ -30,11 +30,7 @@ #include "absl/log/check.h" #include "absl/random/bit_gen_ref.h" #include "absl/types/span.h" -#include "ortools/base/hash.h" -#include "ortools/base/logging.h" -#include "ortools/base/macros.h" #include "ortools/base/strong_vector.h" -#include "ortools/base/types.h" #include "ortools/sat/drat_proof_handler.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" @@ -67,7 +63,12 @@ class SatClause { // Number of literals in the clause. int size() const { return size_; } - int empty() const { return size_ == 0; } + + // We re-use the size to lazily remove clause and notify that they need to be + // deleted. It is why this is not called empty() to emphasis that fact. Note + // that we never create an initially empty clause, so there is no confusion + // with an infeasible model with an empty clause inside. + int IsRemoved() const { return size_ == 0; } // Allows for range based iteration: for (Literal literal : clause) {}. const Literal* begin() const { return &(literals_[0]); } @@ -108,9 +109,6 @@ class SatClause { // be satisfied by completing the assignment. bool IsSatisfied(const VariablesAssignment& assignment) const; - // Returns true if the clause is attached to a LiteralWatchers. - bool IsAttached() const { return size_ > 0; } - std::string DebugString() const; private: @@ -121,7 +119,7 @@ class SatClause { Literal* literals() { return &(literals_[0]); } // Marks the clause so that the next call to CleanUpWatchers() can identify it - // and actually detach it. We use size_ = 0 for this since the clause will + // and actually remove it. We use size_ = 0 for this since the clause will // never be used afterwards. void Clear() { size_ = 0; } @@ -137,8 +135,6 @@ class SatClause { // This class store the literals inline, and literals_ mark the starts of the // variable length portion. Literal literals_[0]; - - DISALLOW_COPY_AND_ASSIGN(SatClause); }; // Clause information used for the clause database management. Note that only @@ -164,6 +160,11 @@ class BinaryImplicationGraph; class LiteralWatchers : public SatPropagator { public: explicit LiteralWatchers(Model* model); + + // This type is neither copyable nor movable. + LiteralWatchers(const LiteralWatchers&) = delete; + LiteralWatchers& operator=(const LiteralWatchers&) = delete; + ~LiteralWatchers() override; // Must be called before adding clauses referring to such variables. @@ -249,7 +250,7 @@ class LiteralWatchers : public SatPropagator { // by the problem clauses and then the learned one that we keep forever. SatClause* NextClauseToMinimize() { for (; to_minimize_index_ < clauses_.size(); ++to_minimize_index_) { - if (!clauses_[to_minimize_index_]->IsAttached()) continue; + if (clauses_[to_minimize_index_]->IsRemoved()) continue; if (!IsRemovable(clauses_[to_minimize_index_])) { return clauses_[to_minimize_index_++]; } @@ -260,6 +261,20 @@ class LiteralWatchers : public SatPropagator { // Restart the scan in NextClauseToMinimize() from the first problem clause. void ResetToMinimizeIndex() { to_minimize_index_ = 0; } + // Really basic algorithm to return a clause to try to probe. We simply + // loop over the clause that we keep forever, in creation order. This starts + // by the problem clauses and then the learned one that we keep forever. + SatClause* NextClauseToProbe() { + for (; to_probe_index_ < clauses_.size(); ++to_probe_index_) { + if (clauses_[to_probe_index_]->IsRemoved()) continue; + return clauses_[to_probe_index_++]; + } + return nullptr; + } + + // Restart the scan in NextClauseToProbe() from the first problem clause. + void ResetToProbeIndex() { to_probe_index_ = 0; } + // During an inprocessing phase, it is easier to detach all clause first, // then simplify and then reattach them. Note however that during these // two calls, it is not possible to use the solver unit-progation. @@ -299,7 +314,7 @@ class LiteralWatchers : public SatPropagator { // Optimization. An index in the clause. Instead of looking for another // literal to watch from the start, we will start from here instead, and - // loop around if needed. This allows to avoid bad quadratric corner cases + // loop around if needed. This allows to avoid bad quadratic corner cases // and lead to an "optimal" complexity. See "Optimal Implementation of // Watched Literals and more General Techniques", Ian P. Gent. // @@ -364,13 +379,12 @@ class LiteralWatchers : public SatPropagator { std::vector clauses_; int to_minimize_index_ = 0; + int to_probe_index_ = 0; // Only contains removable clause. absl::flat_hash_map clauses_info_; DratProofHandler* drat_proof_handler_ = nullptr; - - DISALLOW_COPY_AND_ASSIGN(LiteralWatchers); }; // A binary clause. This is used by BinaryClauseManager. @@ -386,6 +400,11 @@ struct BinaryClause { class BinaryClauseManager { public: BinaryClauseManager() = default; + + // This type is neither copyable nor movable. + BinaryClauseManager(const BinaryClauseManager&) = delete; + BinaryClauseManager& operator=(const BinaryClauseManager&) = delete; + int NumClauses() const { return set_.size(); } // Adds a new binary clause to the manager and returns true if it wasn't @@ -408,7 +427,6 @@ class BinaryClauseManager { private: absl::flat_hash_set> set_; std::vector newly_added_; - DISALLOW_COPY_AND_ASSIGN(BinaryClauseManager); }; // Special class to store and propagate clauses of size 2 (i.e. implication). @@ -462,10 +480,17 @@ class BinaryImplicationGraph : public SatPropagator { stats_("BinaryImplicationGraph"), time_limit_(model->GetOrCreate()), random_(model->GetOrCreate()), - trail_(model->GetOrCreate()) { + trail_(model->GetOrCreate()), + at_most_one_max_expansion_size_( + model->GetOrCreate() + ->at_most_one_max_expansion_size()) { trail_->RegisterPropagator(this); } + // This type is neither copyable nor movable. + BinaryImplicationGraph(const BinaryImplicationGraph&) = delete; + BinaryImplicationGraph& operator=(const BinaryImplicationGraph&) = delete; + ~BinaryImplicationGraph() override { IF_STATS_ENABLED({ LOG(INFO) << stats_.StatString(); @@ -511,9 +536,9 @@ class BinaryImplicationGraph : public SatPropagator { // of literals that will be false if one of the literal in the amo is at one. // It is a way to merge common list of implications. // - // If the final AMO size is smaller than "expansion_size" we fully expand it. - ABSL_MUST_USE_RESULT bool AddAtMostOne(absl::Span at_most_one, - int expansion_size = 10); + // If the final AMO size is smaller than at_most_one_expansion_size + // parameters, we fully expand it. + ABSL_MUST_USE_RESULT bool AddAtMostOne(absl::Span at_most_one); // Uses the binary implication graph to minimize the given conflict by // removing literals that implies others. The idea is that if a and b are two @@ -725,7 +750,9 @@ class BinaryImplicationGraph : public SatPropagator { void RemoveBooleanVariable( BooleanVariable var, std::deque>* postsolve_clauses); bool IsRemoved(Literal l) const { return is_removed_[l]; } - void CleanupAllRemovedVariables(); + void RemoveAllRedundantVariables( + std::deque>* postsolve_clauses); + void CleanupAllRemovedAndFixedVariables(); // ExpandAtMostOneWithWeight() will increase this, so a client can put a limit // on this possibly expansive operation. @@ -739,9 +766,15 @@ class BinaryImplicationGraph : public SatPropagator { const absl::StrongVector& can_be_included, const absl::StrongVector& expanded_lp_values); + // Restarts the at_most_one iterator. + void ResetAtMostOneIterator() { at_most_one_iterator_ = 0; } + + // Returns the next at_most_one, or a span of size 0 when finished. + absl::Span NextAtMostOne(); + private: // Simple wrapper to not forget to output newly fixed variable to the DRAT - // proof if needed. This will propagate rigth away the implications. + // proof if needed. This will propagate right away the implications. bool FixLiteral(Literal true_literal); // Propagates all the direct implications of the given literal becoming true. @@ -769,8 +802,9 @@ class BinaryImplicationGraph : public SatPropagator { // fixed literals and deal with duplicates. Return false iff the model is // UNSAT. // - // If the final AMO size is smaller than "expansion_size" we fully expand it. - bool CleanUpAndAddAtMostOnes(int base_index, int expansion_size = 10); + // If the final AMO size is smaller than the at_most_one_expansion_size + // parameters, we fully expand it. + bool CleanUpAndAddAtMostOnes(int base_index); // To be used in DCHECKs(). bool InvariantsAreOk(); @@ -809,6 +843,8 @@ class BinaryImplicationGraph : public SatPropagator { absl::StrongVector> at_most_ones_; std::vector at_most_one_buffer_; + const int at_most_one_max_expansion_size_; + int at_most_one_iterator_ = 0; // Used by GenerateAtMostOnesWithLargeWeight(). std::vector> tmp_cuts_; @@ -856,8 +892,6 @@ class BinaryImplicationGraph : public SatPropagator { // For RemoveFixedVariables(). int num_processed_fixed_variables_ = 0; - - DISALLOW_COPY_AND_ASSIGN(BinaryImplicationGraph); }; extern template std::vector diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index 7ba5a96064..f089c237a0 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -1903,11 +1903,17 @@ bool IsVarEqOrNeqValue(PresolveContext* context, const LinearConstraintProto& lin) { if (lin.vars_size() != 1) return false; const Domain rhs = ReadDomainFromProto(lin); + + // This is literal => var == value. if (rhs.IsFixed()) return true; - return rhs.InverseMultiplicationBy(lin.coeffs(0)) - .Complement() - .IntersectionWith(context->DomainOf(lin.vars(0))) - .IsFixed(); + + // Is it literal => var != value ? + const Domain not_implied = + rhs.InverseMultiplicationBy(lin.coeffs(0)) + .Complement() + .IntersectionWith(context->DomainOf(lin.vars(0))); + if (not_implied.IsEmpty()) return false; + return not_implied.IsFixed(); } // This method will scan all constraints of all variables appearing in an @@ -2067,6 +2073,7 @@ void MaybeExpandAllDiff(ConstraintProto* ct, PresolveContext* context, context->params().expand_alldiff_constraints(); AllDifferentConstraintProto& proto = *ct->mutable_all_diff(); if (proto.exprs_size() <= 1) return; + if (context->ModelIsUnsat()) return; bool keep_after_expansion = false; bool expand_all_diff_from_usage = false; @@ -2161,6 +2168,7 @@ void ExpandCpModel(PresolveContext* context) { // Make sure all domains are initialized. context->InitializeNewDomains(); + if (context->ModelIsUnsat()) return; // Clear the precedence cache. context->ClearPrecedenceCache(); diff --git a/ortools/sat/cp_model_objective.cc b/ortools/sat/cp_model_objective.cc deleted file mode 100644 index 63151a3b30..0000000000 --- a/ortools/sat/cp_model_objective.cc +++ /dev/null @@ -1,103 +0,0 @@ -// Copyright 2010-2022 Google LLC -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -#include "ortools/sat/cp_model_objective.h" - -#include -#include -#include - -#include "absl/log/check.h" -#include "ortools/sat/cp_model.pb.h" -#include "ortools/sat/cp_model_utils.h" -#include "ortools/util/sorted_interval_list.h" - -namespace operations_research { -namespace sat { - -void EncodeObjectiveAsSingleVariable(CpModelProto* cp_model) { - if (!cp_model->has_objective()) return; - - if (cp_model->objective().vars_size() == 1) { - // Canonicalize the objective to make it easier on us by always making the - // coefficient equal to 1.0. - const int old_ref = cp_model->objective().vars(0); - const int64_t old_coeff = cp_model->objective().coeffs(0); - const double muliplier = static_cast(std::abs(old_coeff)); - if (old_coeff < 0) { - cp_model->mutable_objective()->set_vars(0, NegatedRef(old_ref)); - } - if (muliplier != 1.0) { - // TODO(user): deal with this case. - CHECK(cp_model->objective().domain().empty()); - - double old_factor = cp_model->objective().scaling_factor(); - if (old_factor == 0.0) old_factor = 1.0; - const double old_offset = cp_model->objective().offset(); - cp_model->mutable_objective()->set_scaling_factor(old_factor * muliplier); - cp_model->mutable_objective()->set_offset(old_offset / muliplier); - } - cp_model->mutable_objective()->set_coeffs(0, 1.0); - return; - } - - // Compute trivial bounds on the objective, this is needed otherwise the - // overflow checker might not be happy with the new constraint we are about - // to create. Note that the model validator should make sure that there is no - // overflow in the computation below. - int64_t min_obj = 0; - int64_t max_obj = 0; - for (int i = 0; i < cp_model->objective().vars_size(); ++i) { - const int ref = cp_model->objective().vars(i); - const int var = PositiveRef(ref); - const int64_t coeff = - cp_model->objective().coeffs(i) * (RefIsPositive(ref) ? 1 : -1); - const int64_t value1 = cp_model->variables(var).domain(0) * coeff; - const int64_t value2 = cp_model->variables(var).domain( - cp_model->variables(var).domain_size() - 1) * - coeff; - min_obj += std::min(value1, value2); - max_obj += std::max(value1, value2); - } - - // Create the new objective var. - const int obj_ref = cp_model->variables_size(); - { - IntegerVariableProto* obj = cp_model->add_variables(); - Domain obj_domain(min_obj, max_obj); - if (!cp_model->objective().domain().empty()) { - obj_domain = obj_domain.IntersectionWith( - ReadDomainFromProto(cp_model->objective())); - } - FillDomainInProto(obj_domain, obj); - } - - // Add the linear constraint. - LinearConstraintProto* ct = cp_model->add_constraints()->mutable_linear(); - ct->add_domain(0); - ct->add_domain(0); - *(ct->mutable_vars()) = cp_model->objective().vars(); - *(ct->mutable_coeffs()) = cp_model->objective().coeffs(); - ct->add_vars(obj_ref); - ct->add_coeffs(-1); - - // Update the objective. - cp_model->mutable_objective()->clear_vars(); - cp_model->mutable_objective()->clear_coeffs(); - cp_model->mutable_objective()->add_vars(obj_ref); - cp_model->mutable_objective()->add_coeffs(1); - cp_model->mutable_objective()->clear_domain(); -} - -} // namespace sat -} // namespace operations_research diff --git a/ortools/sat/cp_model_objective.h b/ortools/sat/cp_model_objective.h deleted file mode 100644 index 01928e88fe..0000000000 --- a/ortools/sat/cp_model_objective.h +++ /dev/null @@ -1,40 +0,0 @@ -// Copyright 2010-2022 Google LLC -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -#ifndef OR_TOOLS_SAT_CP_MODEL_OBJECTIVE_H_ -#define OR_TOOLS_SAT_CP_MODEL_OBJECTIVE_H_ - -#include "ortools/sat/cp_model.pb.h" - -namespace operations_research { -namespace sat { - -// If the objective contains more than one term, a new integer variable will be -// created and an equality constraint will be added to the model to force -// this new variable to be equal to the old linear combinaison defining the -// objective. -// -// If the objective contains only one term, we canonicalize it so that its -// coefficient is one. If there is no objecitve, this has no effect. -// -// The main raison to do that is that the objective can be preprocessed like -// any linear constraint before beeing "extracted" again. -// -// TODO(user): Add only a >= constraint? Not sure it is useful to do that before -// presolve, we can always do this optimization just before solving the model. -void EncodeObjectiveAsSingleVariable(CpModelProto* cp_model); - -} // namespace sat -} // namespace operations_research - -#endif // OR_TOOLS_SAT_CP_MODEL_OBJECTIVE_H_ diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index a54f79a606..88a6a3cd8b 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -10560,6 +10560,7 @@ bool CpModelPresolver::ProcessChangedVariables(std::vector* in_queue, for (int i = 0; i < vector_that_can_grow_during_iter.size(); ++i) { const int v = vector_that_can_grow_during_iter[i]; if (context_->VariableIsNotUsedAnymore(v)) continue; + if (context_->ModelIsUnsat()) return false; if (!PresolveAffineRelationIfAny(v)) return false; if (context_->VariableIsNotUsedAnymore(v)) continue; @@ -10688,6 +10689,7 @@ void CpModelPresolver::PresolveToFixPoint() { // Make sure all affine relations are propagated. // This also remove the relation if the degree is now one. + if (context_->ModelIsUnsat()) return; if (!PresolveAffineRelationIfAny(v)) return; const int degree = context_->VarToConstraints(v).size(); diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index 2a6d9eefdb..c953faa9d2 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -605,6 +605,7 @@ absl::flat_hash_map GetNamedParameters( SatParameters new_params = base_params; new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH); new_params.set_use_probing_search(true); + new_params.set_at_most_one_max_expansion_size(2); if (base_params.use_dual_scheduling_heuristics()) { AddDualSchedulingHeuristics(new_params); } diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index be3e19fe7c..11a5be6217 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -2270,7 +2270,7 @@ class FullProblemSolver : public SubSolver { shared_(shared), split_in_chunks_(split_in_chunks), stop_at_first_solution_(stop_at_first_solution), - local_model_(name) { + local_model_(SubSolver::name()) { // Setup the local model parameters and time limit. *(local_model_.GetOrCreate()) = local_parameters; shared_->time_limit->UpdateLocalLimit( diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index 0203e5f27a..c1702ad8c7 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -19,12 +19,15 @@ #include #include #include +#include +#include #include #include #include "absl/container/flat_hash_set.h" #include "absl/log/check.h" #include "absl/types/span.h" +#include "ortools/base/logging.h" #include "ortools/sat/cumulative_energy.h" #include "ortools/sat/diffn_util.h" #include "ortools/sat/disjunctive.h" @@ -191,11 +194,270 @@ void AddNonOverlappingRectangles(const std::vector& x, AddDiffnCumulativeRelationOnX(y_helper, x_helper, model); } } + + if (params.use_area_energetic_reasoning_in_no_overlap_2d()) { + NonOverlappingRectanglesEnergyPropagator* energy_constraint = + new NonOverlappingRectanglesEnergyPropagator(x_helper, y_helper, model); + GenericLiteralWatcher* const watcher = + model->GetOrCreate(); + watcher->SetPropagatorPriority(energy_constraint->RegisterWith(watcher), 5); + model->TakeOwnership(energy_constraint); + } } #define RETURN_IF_FALSE(f) \ if (!(f)) return false; +NonOverlappingRectanglesEnergyPropagator:: + ~NonOverlappingRectanglesEnergyPropagator() { + if (!VLOG_IS_ON(1)) return; + std::vector> stats; + stats.push_back( + {"NonOverlappingRectanglesEnergyPropagator/called", num_calls_}); + stats.push_back( + {"NonOverlappingRectanglesEnergyPropagator/conflicts", num_conflicts_}); + stats.push_back( + {"NonOverlappingRectanglesEnergyPropagator/multiple_conflicts", + num_multiple_conflicts_}); + shared_stats_->AddStats(stats); +} + +bool NonOverlappingRectanglesEnergyPropagator::Propagate() { + // TODO(user): double-check/revisit the algo for box of variable sizes. + const int num_boxes = x_.NumTasks(); + if (!x_.SynchronizeAndSetTimeDirection(true)) return false; + if (!y_.SynchronizeAndSetTimeDirection(true)) return false; + + num_calls_++; + Rectangle bounding_box = {.x_min = std::numeric_limits::max(), + .x_max = std::numeric_limits::min(), + .y_min = std::numeric_limits::max(), + .y_max = std::numeric_limits::min()}; + std::vector active_box_ranges; + active_box_ranges.reserve(num_boxes); + for (int box = 0; box < num_boxes; ++box) { + if (x_.SizeMin(box) == 0 || y_.SizeMin(box) == 0) continue; + if (!x_.IsPresent(box) || !y_.IsPresent(box)) continue; + + bounding_box.x_min = std::min(bounding_box.x_min, x_.StartMin(box)); + bounding_box.x_max = std::max(bounding_box.x_max, x_.EndMax(box)); + bounding_box.y_min = std::min(bounding_box.y_min, y_.StartMin(box)); + bounding_box.y_max = std::max(bounding_box.y_max, y_.EndMax(box)); + + active_box_ranges.push_back(RectangleInRange{ + .box_index = box, + .bounding_area = {.x_min = x_.StartMin(box), + .x_max = x_.StartMax(box) + x_.SizeMin(box), + .y_min = y_.StartMin(box), + .y_max = y_.StartMax(box) + y_.SizeMin(box)}, + .x_size = x_.SizeMin(box), + .y_size = y_.SizeMin(box)}); + } + + if (active_box_ranges.size() <= 2) { + return true; + } + + // Our algo is quadratic, so we don't want to run it on really large problems. + if (active_box_ranges.size() > 1000) { + return true; + } + + if (std::max(bounding_box.SizeX(), bounding_box.SizeY()) * + active_box_ranges.size() > + std::numeric_limits::max()) { + // Avoid integer overflows if the area of the boxes get comparable with + // INT64_MAX + return true; + } + + const std::vector rectangles_too_much_energy = + FindRectanglesWithEnergyConflictMC(active_box_ranges, *random_, 1.0); + + if (rectangles_too_much_energy.empty()) return true; + + num_conflicts_++; + num_multiple_conflicts_ += rectangles_too_much_energy.size() > 1; + + std::optional> best_explanation; + for (const auto& r : rectangles_too_much_energy) { + std::vector range_for_explanation = + GetEnergyConflictForRectangle(r, active_box_ranges); + CheckPropagationIsValid(range_for_explanation, r); + if (!best_explanation || + best_explanation->size() > range_for_explanation.size()) { + best_explanation = range_for_explanation; + } + } + BuildAndReportEnergyTooLarge(*best_explanation); + return false; +} + +std::vector +NonOverlappingRectanglesEnergyPropagator::GetEnergyConflictForRectangle( + const Rectangle& rectangle, + const std::vector& active_box_ranges) { + struct OverlapPerBox { + IntegerValue energy; + IntegerValue overlap_x_size; + IntegerValue overlap_y_size; + int index; + }; + + // First select the minimum amount of elements that would be enough to exceed + // the energy. + std::vector energy_per_box(active_box_ranges.size()); + for (int i = 0; i < active_box_ranges.size(); ++i) { + const Rectangle overlap = + active_box_ranges[i].GetMinimumIntersection(rectangle); + OverlapPerBox& box = energy_per_box[i]; + box.overlap_x_size = overlap.SizeX(); + box.overlap_y_size = overlap.SizeY(); + box.index = i; + box.energy = box.overlap_x_size * box.overlap_y_size; + } + std::sort(energy_per_box.begin(), energy_per_box.end(), + [](const OverlapPerBox& a, const OverlapPerBox& b) { + return a.energy > b.energy; + }); + IntegerValue available_energy = rectangle.Area(); + IntegerValue used_energy = 0; + std::vector ranges_for_explanation; + ranges_for_explanation.reserve(energy_per_box.size()); + for (int i = 0; i < energy_per_box.size(); ++i) { + const OverlapPerBox& box = energy_per_box[i]; + used_energy += box.energy; + ranges_for_explanation.push_back( + RectangleInRange::BiggestWithMinIntersection( + rectangle, active_box_ranges[box.index], box.overlap_x_size, + box.overlap_y_size)); + if (used_energy > available_energy) { + break; + } + } + + // Now use the energy slack to make the intervals even bigger. + IntegerValue slack = used_energy - available_energy - 1; + VLOG_EVERY_N_SEC(2, 3) << "Rectangle with energy overflow: " << rectangle + << " (" << used_energy << "/" << available_energy + << "), with " << ranges_for_explanation.size() + << " boxes inside."; + + for (auto& range : ranges_for_explanation) { + Rectangle overlap = range.GetMinimumIntersection(rectangle); + DCHECK_GT(overlap.SizeX(), 0); + DCHECK_GT(overlap.SizeY(), 0); + if (overlap.SizeX() < slack) { + IntegerValue y_slack = slack / overlap.SizeX(); + // We know we can remove an area as big as y_slack * x_size. That means + // we can reduce the start of an interval of y_slack or increase the end + // by the same amount. But there is no point wasting the slack on a + // interval larger than the zero-level interval, so we check for that. + const IntegerValue y_min_slack = std::max( + IntegerValue(0), + std::min(y_slack, range.bounding_area.y_min - + y_.LevelZeroStartMin(range.box_index))); + if (y_min_slack > 0) { + range.bounding_area.y_min -= y_min_slack; + slack -= y_min_slack * overlap.SizeX(); + y_slack -= y_min_slack; + } + const IntegerValue y_max_slack = + std::min(y_slack, y_.LevelZeroStartMax(range.box_index) - + range.bounding_area.y_max); + if (y_max_slack > 0) { + range.bounding_area.y_max += y_max_slack; + slack -= y_max_slack * overlap.SizeX(); + } + // Recompute the overlap, since we changed the item dimensions. + overlap = range.GetMinimumIntersection(rectangle); + } + if (overlap.SizeY() < slack) { + IntegerValue x_slack = slack / overlap.SizeY(); + const IntegerValue x_min_slack = std::max( + IntegerValue(0), + std::min(x_slack, range.bounding_area.x_min - + x_.LevelZeroStartMax(range.box_index))); + if (x_min_slack > 0) { + range.bounding_area.x_min -= x_min_slack; + slack -= x_min_slack * overlap.SizeY(); + x_slack -= x_min_slack; + } + const IntegerValue x_max_slack = + std::min(x_slack, x_.LevelZeroStartMax(range.box_index) - + range.bounding_area.x_max); + if (x_max_slack > 0) { + range.bounding_area.x_max += x_max_slack; + slack -= x_max_slack * overlap.SizeY(); + } + } + } + CHECK_GT(used_energy, available_energy); + CHECK_GT(available_energy, 0); + return ranges_for_explanation; +} + +int NonOverlappingRectanglesEnergyPropagator::RegisterWith( + GenericLiteralWatcher* watcher) { + const int id = watcher->Register(this); + x_.WatchAllTasks(id, watcher, /*watch_start_max=*/true, + /*watch_end_max=*/true); + y_.WatchAllTasks(id, watcher, /*watch_start_max=*/true, + /*watch_end_max=*/true); + return id; +} + +void NonOverlappingRectanglesEnergyPropagator::CheckPropagationIsValid( + const std::vector& ranges, + const Rectangle& rectangle_too_much_energy) { + const IntegerValue available_energy = rectangle_too_much_energy.Area(); + IntegerValue used_energy = 0; + for (const auto& range : ranges) { + const int b = range.box_index; + // Check that we are explaining intervals at least as large as the current. + CHECK_GE(x_.StartMin(b), range.bounding_area.x_min); + CHECK_GE(range.bounding_area.x_max - range.x_size, x_.StartMax(b)); + CHECK_GE(y_.StartMin(b), range.bounding_area.y_min); + CHECK_GE(range.bounding_area.y_max - range.y_size, y_.StartMax(b)); + + // Each one of the boxes-in-range that we found on the cut does intersect + // the rectangle we found. + const auto intersection = + range.GetMinimumIntersection(rectangle_too_much_energy); + CHECK_GT(intersection.Area(), 0); + // It cannot intersect more than the size of the object. + CHECK_GE(x_.SizeMin(b), intersection.SizeX()); + CHECK_GE(y_.SizeMin(b), intersection.SizeY()); + used_energy += intersection.Area(); + } + + // We have more area inside the rectangle than the area of the rectangle. + CHECK_GT(used_energy, available_energy); +} + +bool NonOverlappingRectanglesEnergyPropagator::BuildAndReportEnergyTooLarge( + const std::vector& ranges) { + x_.ClearReason(); + y_.ClearReason(); + for (const auto& range : ranges) { + const int b = range.box_index; + + x_.AddStartMinReason(b, range.bounding_area.x_min); + y_.AddStartMinReason(b, range.bounding_area.y_min); + + x_.AddStartMaxReason(b, range.bounding_area.x_max - range.x_size); + y_.AddStartMaxReason(b, range.bounding_area.y_max - range.y_size); + + x_.AddSizeMinReason(b); + y_.AddSizeMinReason(b); + + x_.AddPresenceReason(b); + y_.AddPresenceReason(b); + } + x_.ImportOtherReasons(y_); + return x_.ReportConflict(); +} + namespace { // We want for different propagation to reuse as much as possible the same diff --git a/ortools/sat/diffn.h b/ortools/sat/diffn.h index 9242f2499e..a4e3f2c0db 100644 --- a/ortools/sat/diffn.h +++ b/ortools/sat/diffn.h @@ -14,6 +14,7 @@ #ifndef OR_TOOLS_SAT_DIFFN_H_ #define OR_TOOLS_SAT_DIFFN_H_ +#include #include #include @@ -26,10 +27,52 @@ #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" #include "ortools/sat/sat_parameters.pb.h" +#include "ortools/sat/synchronization.h" +#include "ortools/sat/util.h" namespace operations_research { namespace sat { +// Propagates using a box energy reasoning. +class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface { + public: + NonOverlappingRectanglesEnergyPropagator(SchedulingConstraintHelper* x, + SchedulingConstraintHelper* y, + Model* model) + : x_(*x), + y_(*y), + random_(model->GetOrCreate()), + shared_stats_(model->GetOrCreate()) {} + + ~NonOverlappingRectanglesEnergyPropagator() override; + + bool Propagate() final; + int RegisterWith(GenericLiteralWatcher* watcher); + + private: + bool BuildAndReportEnergyTooLarge( + const std::vector& ranges); + void CheckPropagationIsValid(const std::vector& ranges, + const Rectangle& rectangle_too_much_energy); + std::vector GetEnergyConflictForRectangle( + const Rectangle& rectangle, + const std::vector& active_box_ranges); + + SchedulingConstraintHelper& x_; + SchedulingConstraintHelper& y_; + ModelRandomGenerator* random_; + SharedStatistics* shared_stats_; + + int64_t num_calls_ = 0; + int64_t num_conflicts_ = 0; + int64_t num_multiple_conflicts_ = 0; + + NonOverlappingRectanglesEnergyPropagator( + const NonOverlappingRectanglesEnergyPropagator&) = delete; + NonOverlappingRectanglesEnergyPropagator& operator=( + const NonOverlappingRectanglesEnergyPropagator&) = delete; +}; + // Enforces that the boxes with corners in (x, y), (x + dx, y), (x, y + dy) // and (x + dx, y + dy) do not overlap. void AddNonOverlappingRectangles(const std::vector& x, diff --git a/ortools/sat/diffn_util.cc b/ortools/sat/diffn_util.cc index 2f881b0671..cbfba9976d 100644 --- a/ortools/sat/diffn_util.cc +++ b/ortools/sat/diffn_util.cc @@ -16,13 +16,16 @@ #include #include +#include #include #include #include #include "absl/container/flat_hash_set.h" +#include "absl/container/inlined_vector.h" #include "absl/log/check.h" #include "absl/random/bit_gen_ref.h" +#include "absl/random/discrete_distribution.h" #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" @@ -604,5 +607,403 @@ IntegerValue CapacityProfile::GetBoundingArea() { return area; } +IntegerValue Smallest1DIntersection(IntegerValue range_min, + IntegerValue range_max, IntegerValue size, + IntegerValue interval_min, + IntegerValue interval_max) { + // If the item is on the left of the range, we get the intersection between + // [range_min, range_min + size] and [interval_min, interval_max]. + const IntegerValue overlap_on_left = + std::min(range_min + size, interval_max) - + std::max(range_min, interval_min); + + // If the item is on the right of the range, we get the intersection between + // [range_max - size, range_max] and [interval_min, interval_max]. + const IntegerValue overlap_on_right = + std::min(range_max, interval_max) - + std::max(range_max - size, interval_min); + + return std::max(IntegerValue(0), std::min(overlap_on_left, overlap_on_right)); +} + +ProbingRectangle::ProbingRectangle( + const std::vector& intervals) + : intervals_(intervals) { + minimum_energy_ = 0; + if (intervals_.empty()) { + return; + } + interval_points_sorted_by_x_.reserve(intervals_.size() * 4); + interval_points_sorted_by_y_.reserve(intervals_.size() * 4); + for (int i = 0; i < intervals_.size(); ++i) { + const RectangleInRange& interval = intervals_[i]; + minimum_energy_ += interval.x_size * interval.y_size; + + interval_points_sorted_by_x_.push_back( + {interval.bounding_area.x_min, i, + IntervalPoint::IntervalPointType::START_MIN}); + interval_points_sorted_by_x_.push_back( + {interval.bounding_area.x_min + interval.x_size, i, + IntervalPoint::IntervalPointType::END_MIN}); + interval_points_sorted_by_x_.push_back( + {interval.bounding_area.x_max - interval.x_size, i, + IntervalPoint::IntervalPointType::START_MAX}); + interval_points_sorted_by_x_.push_back( + {interval.bounding_area.x_max, i, + IntervalPoint::IntervalPointType::END_MAX}); + + interval_points_sorted_by_y_.push_back( + {interval.bounding_area.y_min, i, + IntervalPoint::IntervalPointType::START_MIN}); + interval_points_sorted_by_y_.push_back( + {interval.bounding_area.y_min + interval.y_size, i, + IntervalPoint::IntervalPointType::END_MIN}); + interval_points_sorted_by_y_.push_back( + {interval.bounding_area.y_max - interval.y_size, i, + IntervalPoint::IntervalPointType::START_MAX}); + interval_points_sorted_by_y_.push_back( + {interval.bounding_area.y_max, i, + IntervalPoint::IntervalPointType::END_MAX}); + } + + std::sort(interval_points_sorted_by_x_.begin(), + interval_points_sorted_by_x_.end(), + [](const IntervalPoint& a, const IntervalPoint& b) { + return a.value < b.value; + }); + std::sort(interval_points_sorted_by_y_.begin(), + interval_points_sorted_by_y_.end(), + [](const IntervalPoint& a, const IntervalPoint& b) { + return a.value < b.value; + }); + + grouped_intervals_sorted_by_x_.reserve(interval_points_sorted_by_x_.size()); + grouped_intervals_sorted_by_y_.reserve(interval_points_sorted_by_y_.size()); + int i = 0; + while (i < interval_points_sorted_by_x_.size()) { + int idx_begin = i; + while (i < interval_points_sorted_by_x_.size() && + interval_points_sorted_by_x_[i].value == + interval_points_sorted_by_x_[idx_begin].value) { + i++; + } + grouped_intervals_sorted_by_x_.push_back( + {interval_points_sorted_by_x_[idx_begin].value, + absl::Span(interval_points_sorted_by_x_) + .subspan(idx_begin, i - idx_begin)}); + } + + i = 0; + while (i < interval_points_sorted_by_y_.size()) { + int idx_begin = i; + while (i < interval_points_sorted_by_y_.size() && + interval_points_sorted_by_y_[i].value == + interval_points_sorted_by_y_[idx_begin].value) { + i++; + } + grouped_intervals_sorted_by_y_.push_back( + {interval_points_sorted_by_y_[idx_begin].value, + absl::Span(interval_points_sorted_by_y_) + .subspan(idx_begin, i - idx_begin)}); + } + + left_index_ = 0; + right_index_ = grouped_intervals_sorted_by_x_.size() - 1; + bottom_index_ = 0; + top_index_ = grouped_intervals_sorted_by_y_.size() - 1; + + for (const auto& point : grouped_intervals_sorted_by_x_[left_index_].points) { + ranges_touching_boundary_[Edge::LEFT].insert(point.index); + } + for (const auto& point : + grouped_intervals_sorted_by_x_[right_index_].points) { + ranges_touching_boundary_[Edge::RIGHT].insert(point.index); + } + for (const auto& point : + grouped_intervals_sorted_by_y_[bottom_index_].points) { + ranges_touching_boundary_[Edge::BOTTOM].insert(point.index); + } + for (const auto& point : grouped_intervals_sorted_by_y_[top_index_].points) { + ranges_touching_boundary_[Edge::TOP].insert(point.index); + } + probe_area_ = GetCurrentRectangle().Area(); +} + +Rectangle ProbingRectangle::GetCurrentRectangle() const { + return {.x_min = grouped_intervals_sorted_by_x_[left_index_].coordinate, + .x_max = grouped_intervals_sorted_by_x_[right_index_].coordinate, + .y_min = grouped_intervals_sorted_by_y_[bottom_index_].coordinate, + .y_max = grouped_intervals_sorted_by_y_[top_index_].coordinate}; +} + +void ProbingRectangle::Shrink(Edge edge) { + absl::Span points; + + minimum_energy_ -= GetShrinkDeltaEnergy(edge); + switch (edge) { + case Edge::LEFT: + left_index_++; + points = grouped_intervals_sorted_by_x_[left_index_].points; + break; + case Edge::BOTTOM: + bottom_index_++; + points = grouped_intervals_sorted_by_y_[bottom_index_].points; + break; + case Edge::RIGHT: + right_index_--; + points = grouped_intervals_sorted_by_x_[right_index_].points; + break; + case Edge::TOP: + top_index_--; + points = grouped_intervals_sorted_by_y_[top_index_].points; + break; + } + + for (const auto& point : points) { + const bool became_outside_probe = + (point.type == IntervalPoint::IntervalPointType::END_MIN && + (edge == Edge::LEFT || edge == Edge::BOTTOM)) || + (point.type == IntervalPoint::IntervalPointType::START_MAX && + (edge == Edge::RIGHT || edge == Edge::TOP)); + if (became_outside_probe) { + ranges_touching_boundary_[Edge::LEFT].erase(point.index); + ranges_touching_boundary_[Edge::BOTTOM].erase(point.index); + ranges_touching_boundary_[Edge::RIGHT].erase(point.index); + ranges_touching_boundary_[Edge::TOP].erase(point.index); + } + } + + const Rectangle current_rectangle = GetCurrentRectangle(); + auto can_consume_energy = [¤t_rectangle]( + const RectangleInRange& range) { + // This intersects the current rectangle with the largest rectangle + // that must intersect with the range in some way. To visualize this + // largest rectangle, imagine the four possible extreme positions for + // the item in range (the four corners). This rectangle is the one + // defined by the interior points of each position. + // This don't use IsDisjoint() because it also works when the rectangle + // would be malformed (it's bounding box less than twice the size). + return !( + range.bounding_area.x_max - range.x_size >= current_rectangle.x_max || + range.bounding_area.y_max - range.y_size >= current_rectangle.y_max || + current_rectangle.x_min >= range.bounding_area.x_min + range.x_size || + current_rectangle.y_min >= range.bounding_area.y_min + range.y_size); + }; + + switch (edge) { + case Edge::LEFT: + case Edge::BOTTOM: + for (const auto& point : points) { + if (point.type == IntervalPoint::IntervalPointType::START_MIN) { + if (can_consume_energy(intervals_[point.index])) { + ranges_touching_boundary_[edge].insert(point.index); + } + } + } + break; + case Edge::RIGHT: + case Edge::TOP: + for (const auto& point : points) { + if (point.type == IntervalPoint::IntervalPointType::END_MAX) { + if (can_consume_energy(intervals_[point.index])) { + ranges_touching_boundary_[edge].insert(point.index); + } + } + } + break; + } + probe_area_ = GetCurrentRectangle().Area(); +} + +IntegerValue ProbingRectangle::GetShrinkDeltaArea(Edge edge) const { + const Rectangle current_rectangle = GetCurrentRectangle(); + switch (edge) { + case Edge::LEFT: + return (grouped_intervals_sorted_by_x_[left_index_ + 1].coordinate - + current_rectangle.x_min) * + current_rectangle.SizeY(); + case Edge::BOTTOM: + return (grouped_intervals_sorted_by_y_[bottom_index_ + 1].coordinate - + current_rectangle.y_min) * + current_rectangle.SizeX(); + case Edge::RIGHT: + return (current_rectangle.x_max - + grouped_intervals_sorted_by_x_[right_index_ - 1].coordinate) * + current_rectangle.SizeY(); + case Edge::TOP: + return (current_rectangle.y_max - + grouped_intervals_sorted_by_y_[top_index_ - 1].coordinate) * + current_rectangle.SizeX(); + } +} + +IntegerValue ProbingRectangle::GetShrinkDeltaEnergy(Edge edge) const { + const Rectangle current_rectangle = GetCurrentRectangle(); + Rectangle next_rectangle = current_rectangle; + IntegerValue step_1d_size; + + switch (edge) { + case Edge::LEFT: + next_rectangle.x_min = + grouped_intervals_sorted_by_x_[left_index_ + 1].coordinate; + step_1d_size = next_rectangle.x_min - current_rectangle.x_min; + break; + case Edge::BOTTOM: + next_rectangle.y_min = + grouped_intervals_sorted_by_y_[bottom_index_ + 1].coordinate; + step_1d_size = next_rectangle.y_min - current_rectangle.y_min; + break; + case Edge::RIGHT: + next_rectangle.x_max = + grouped_intervals_sorted_by_x_[right_index_ - 1].coordinate; + step_1d_size = current_rectangle.x_max - next_rectangle.x_max; + break; + case Edge::TOP: + next_rectangle.y_max = + grouped_intervals_sorted_by_y_[top_index_ - 1].coordinate; + step_1d_size = current_rectangle.y_max - next_rectangle.y_max; + break; + } + + IntegerValue delta_energy = 0; + IntegerValue units_crossed = 0; + // Note that the non-deterministic iteration order is fine here. + for (const int idx : ranges_touching_boundary_[edge]) { + const RectangleInRange& range = intervals_[idx]; + bool problematic_case_in_two_sides = false; + IntegerValue opposite_slack; + switch (edge) { + case Edge::LEFT: + opposite_slack = range.bounding_area.x_max - current_rectangle.x_max; + // First check if we touch the opposite edge to the one we are + // shrinking. + if (opposite_slack >= 0) { + // If it do, it's problematic if it has more slack on the opposite + // side so it will "jump" to the other side. + problematic_case_in_two_sides = + opposite_slack >= + current_rectangle.x_min - range.bounding_area.x_min; + } + break; + case Edge::BOTTOM: + opposite_slack = range.bounding_area.y_max - current_rectangle.y_max; + if (opposite_slack >= 0) { + problematic_case_in_two_sides = + opposite_slack >= + current_rectangle.y_min - range.bounding_area.y_min; + } + break; + case Edge::RIGHT: + opposite_slack = current_rectangle.x_min - range.bounding_area.x_min; + if (opposite_slack >= 0) { + problematic_case_in_two_sides = + opposite_slack >= + range.bounding_area.x_max - current_rectangle.x_max; + } + break; + case Edge::TOP: + opposite_slack = current_rectangle.y_min - range.bounding_area.y_min; + if (opposite_slack >= 0) { + problematic_case_in_two_sides = + opposite_slack >= + range.bounding_area.y_max - current_rectangle.y_max; + } + break; + } + if (problematic_case_in_two_sides) { + // When it touches both sides, reducing the probe on the bottom might + // make the place with the minimum overlap become the top. It's too + // complicated to manage, so we fall back on actually computing it from + // scratch. + delta_energy += range.GetMinimumIntersectionArea(current_rectangle); + delta_energy -= range.GetMinimumIntersectionArea(next_rectangle); + } else { + IntegerValue intersect_length; + if (edge == Edge::LEFT || edge == Edge::RIGHT) { + intersect_length = Smallest1DIntersection( + range.bounding_area.y_min, range.bounding_area.y_max, range.y_size, + current_rectangle.y_min, current_rectangle.y_max); + } else { + intersect_length = Smallest1DIntersection( + range.bounding_area.x_min, range.bounding_area.x_max, range.x_size, + current_rectangle.x_min, current_rectangle.x_max); + } + units_crossed += intersect_length; + } + } + delta_energy += units_crossed * step_1d_size; + return delta_energy; +} + +bool ProbingRectangle::CanShrink(Edge edge) const { + switch (edge) { + case Edge::LEFT: + case Edge::RIGHT: + return (right_index_ - left_index_ > 1); + case Edge::BOTTOM: + case Edge::TOP: + return (top_index_ - bottom_index_ > 1); + } +} + +namespace { +std::vector GetExpTable() { + std::vector table(101); + for (int i = 0; i <= 100; ++i) { + table[i] = std::exp(-(i - 50) / 5.0); + } + return table; +} +} // namespace + +std::vector FindRectanglesWithEnergyConflictMC( + const std::vector& intervals, absl::BitGenRef random, + double temperature) { + std::vector result; + ProbingRectangle ranges(intervals); + + static const std::vector* cached_probabilities = + new std::vector(GetExpTable()); + + const double inv_temp = 1.0 / temperature; + absl::InlinedVector candidates; + absl::InlinedVector weights; + while (!ranges.IsMinimal()) { + const IntegerValue rect_area = ranges.GetCurrentRectangleArea(); + const IntegerValue min_energy = ranges.GetMinimumEnergy(); + if (min_energy > rect_area) { + result.push_back(ranges.GetCurrentRectangle()); + } + if (min_energy == 0) { + break; + } + candidates.clear(); + weights.clear(); + + for (int border_idx = 0; border_idx < 4; ++border_idx) { + const ProbingRectangle::Edge border = + static_cast(border_idx); + if (!ranges.CanShrink(border)) { + continue; + } + candidates.push_back(border); + const IntegerValue delta_area = ranges.GetShrinkDeltaArea(border); + const IntegerValue delta_energy = ranges.GetShrinkDeltaEnergy(border); + const IntegerValue delta_slack = delta_energy - delta_area; + const int table_lookup = std::max( + 0, std::min((int)(delta_slack.value() * 5 * inv_temp + 50), 100)); + weights.push_back(cached_probabilities->at(table_lookup)); + } + // Pick a change with a probability proportional to exp(- delta_E / Temp) + absl::discrete_distribution dist(weights.begin(), weights.end()); + ranges.Shrink(candidates.at(dist(random))); + } + CHECK_GT(ranges.GetCurrentRectangleArea(), 0); + if (ranges.GetMinimumEnergy() > ranges.GetCurrentRectangleArea()) { + result.push_back(ranges.GetCurrentRectangle()); + } + return result; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/diffn_util.h b/ortools/sat/diffn_util.h index ad0bfba408..6a35661385 100644 --- a/ortools/sat/diffn_util.h +++ b/ortools/sat/diffn_util.h @@ -15,18 +15,19 @@ #define OR_TOOLS_SAT_DIFFN_UTIL_H_ #include -#include #include +#include +#include #include -#include #include +#include #include #include "absl/container/flat_hash_set.h" +#include "absl/log/check.h" #include "absl/random/bit_gen_ref.h" #include "absl/strings/str_format.h" #include "absl/types/span.h" -#include "ortools/graph/connected_components.h" #include "ortools/sat/integer.h" #include "ortools/sat/intervals.h" #include "ortools/util/strong_integers.h" @@ -47,16 +48,65 @@ struct Rectangle { y_max = std::max(y_max, other.y_max); } - IntegerValue Area() const { return (x_max - x_min) * (y_max - y_min); } + IntegerValue Area() const { return SizeX() * SizeY(); } + + IntegerValue SizeX() const { return x_max - x_min; } + IntegerValue SizeY() const { return y_max - y_min; } bool IsDisjoint(const Rectangle& other) const; - std::string DebugString() const { - return absl::StrFormat("rectangle(x(%i..%i), y(%i..%i))", x_min.value(), - x_max.value(), y_min.value(), y_max.value()); + // Returns an empty rectangle if no intersection. + Rectangle Intersect(const Rectangle& other) const; + IntegerValue IntersectArea(const Rectangle& other) const; + + template + friend void AbslStringify(Sink& sink, const Rectangle& r) { + absl::Format(&sink, "rectangle(x(%i..%i), y(%i..%i))", r.x_min.value(), + r.x_max.value(), r.y_min.value(), r.y_max.value()); + } + + bool operator==(const Rectangle& other) const { + return std::tie(x_min, x_max, y_min, y_max) == + std::tie(other.x_min, other.x_max, other.y_min, other.y_max); + } + + static Rectangle GetEmpty() { + return Rectangle{.x_min = IntegerValue(0), + .x_max = IntegerValue(0), + .y_min = IntegerValue(0), + .y_max = IntegerValue(0)}; } }; +inline Rectangle Rectangle::Intersect(const Rectangle& other) const { + const IntegerValue ret_x_min = std::max(x_min, other.x_min); + const IntegerValue ret_y_min = std::max(y_min, other.y_min); + const IntegerValue ret_x_max = std::min(x_max, other.x_max); + const IntegerValue ret_y_max = std::min(y_max, other.y_max); + + if (ret_x_min >= ret_x_max || ret_y_min >= ret_y_max) { + return GetEmpty(); + } else { + return Rectangle{.x_min = ret_x_min, + .x_max = ret_x_max, + .y_min = ret_y_min, + .y_max = ret_y_max}; + } +} + +inline IntegerValue Rectangle::IntersectArea(const Rectangle& other) const { + const IntegerValue ret_x_min = std::max(x_min, other.x_min); + const IntegerValue ret_y_min = std::max(y_min, other.y_min); + const IntegerValue ret_x_max = std::min(x_max, other.x_max); + const IntegerValue ret_y_max = std::min(y_max, other.y_max); + + if (ret_x_min >= ret_x_max || ret_y_min >= ret_y_max) { + return 0; + } else { + return (ret_x_max - ret_x_min) * (ret_y_max - ret_y_min); + } +} + // Creates a graph when two nodes are connected iff their rectangles overlap. // Then partition into connected components. // @@ -255,6 +305,196 @@ class CapacityProfile { int num_rectangles_added_ = 0; }; +// 1D counterpart of RectangleInRange::GetMinimumIntersectionArea. +// Finds the minimum possible overlap of a interval of size `size` that fits in +// [range_min, range_max] and a second interval [interval_min, interval_max]. +IntegerValue Smallest1DIntersection(IntegerValue range_min, + IntegerValue range_max, IntegerValue size, + IntegerValue interval_min, + IntegerValue interval_max); + +// A rectangle of size (`x_size`, `y_size`) that can freely move inside the +// `bounding_area` rectangle. +struct RectangleInRange { + int box_index; + Rectangle bounding_area; + IntegerValue x_size; + IntegerValue y_size; + + enum class Corner { + BOTTOM_LEFT = 0, + TOP_LEFT = 1, + BOTTOM_RIGHT = 2, + TOP_RIGHT = 3, + }; + + // Returns the position of the rectangle fixed to one of the corner of its + // range. + Rectangle GetAtCorner(Corner p) const { + switch (p) { + case Corner::BOTTOM_LEFT: + return Rectangle{.x_min = bounding_area.x_min, + .x_max = bounding_area.x_min + x_size, + .y_min = bounding_area.y_min, + .y_max = bounding_area.y_min + y_size}; + case Corner::TOP_LEFT: + return Rectangle{.x_min = bounding_area.x_min, + .x_max = bounding_area.x_min + x_size, + .y_min = bounding_area.y_max - y_size, + .y_max = bounding_area.y_max}; + case Corner::BOTTOM_RIGHT: + return Rectangle{.x_min = bounding_area.x_max - x_size, + .x_max = bounding_area.x_max, + .y_min = bounding_area.y_min, + .y_max = bounding_area.y_min + y_size}; + case Corner::TOP_RIGHT: + return Rectangle{.x_min = bounding_area.x_max - x_size, + .x_max = bounding_area.x_max, + .y_min = bounding_area.y_max - y_size, + .y_max = bounding_area.y_max}; + } + } + + // Returns an empty rectangle if it is possible for no intersection to happen. + Rectangle GetMinimumIntersection(const Rectangle& containing_area) const { + IntegerValue smallest_area = std::numeric_limits::max(); + Rectangle best_intersection; + for (int corner_idx = 0; corner_idx < 4; ++corner_idx) { + const Corner p = static_cast(corner_idx); + Rectangle intersection = containing_area.Intersect(GetAtCorner(p)); + const IntegerValue intersection_area = intersection.Area(); + if (intersection_area == 0) { + return Rectangle::GetEmpty(); + } + if (intersection_area < smallest_area) { + smallest_area = intersection_area; + best_intersection = std::move(intersection); + } + } + return best_intersection; + } + + IntegerValue GetMinimumIntersectionArea( + const Rectangle& containing_area) const { + return Smallest1DIntersection(bounding_area.x_min, bounding_area.x_max, + x_size, containing_area.x_min, + containing_area.x_max) * + Smallest1DIntersection(bounding_area.y_min, bounding_area.y_max, + y_size, containing_area.y_min, + containing_area.y_max); + } + + static RectangleInRange BiggestWithMinIntersection( + const Rectangle& containing_area, const RectangleInRange& original, + const IntegerValue& min_intersect_x, + const IntegerValue& min_intersect_y) { + const IntegerValue x_size = original.x_size; + const IntegerValue y_size = original.y_size; + + RectangleInRange result; + result.x_size = x_size; + result.y_size = y_size; + result.box_index = original.box_index; + + // We cannot intersect more units than the whole item. + DCHECK_GE(x_size, min_intersect_x); + DCHECK_GE(y_size, min_intersect_y); + + // Units that can *not* intersect per dimension. + const IntegerValue x_headroom = x_size - min_intersect_x; + const IntegerValue y_headroom = y_size - min_intersect_y; + + result.bounding_area.x_min = containing_area.x_min - x_headroom; + result.bounding_area.x_max = containing_area.x_max + x_headroom; + result.bounding_area.y_min = containing_area.y_min - y_headroom; + result.bounding_area.y_max = containing_area.y_max + y_headroom; + + return result; + } +}; + +// Cheaply test several rectangles for area conflict. +// This is used by FindRectanglesWithEnergyConflictMC() below. +class ProbingRectangle { + public: + // It will initialize with the bounding box of the whole set. + explicit ProbingRectangle(const std::vector& intervals); + + enum Edge { TOP = 0, LEFT = 1, BOTTOM = 2, RIGHT = 3 }; + + // Shrink the rectangle by moving one of its four edges to the next + // "interesting" value. The interesting values for x or y are the ones that + // correspond to a boundary, ie., a value that corresponds to one of {min, + // min + size, max - size, max} of a rectangle. + void Shrink(Edge edge); + + bool CanShrink(Edge edge) const; + + bool IsMinimal() const { + // We only need to know if there is slack on both dimensions. Actually + // CanShrink(BOTTOM) == CanShrink(TOP) and conversely. + return !(CanShrink(Edge::BOTTOM) || CanShrink(Edge::LEFT)); + } + + // How much of GetMinimumEnergy() will change if Shrink() is called. + IntegerValue GetShrinkDeltaEnergy(Edge edge) const; + + // How much of GetCurrentRectangleArea() will change if Shrink() is called. + IntegerValue GetShrinkDeltaArea(Edge edge) const; + + Rectangle GetCurrentRectangle() const; + IntegerValue GetCurrentRectangleArea() const { return probe_area_; } + + // This is equivalent of, for every item: + // - Call GetMinimumIntersectionArea() with GetCurrentRectangle(). + // - Return the total sum of the areas. + IntegerValue GetMinimumEnergy() const { return minimum_energy_; } + const std::vector& Intervals() const { return intervals_; } + + private: + struct IntervalPoint { + IntegerValue value; + int index; + enum class IntervalPointType { + START_MIN, + START_MAX, + END_MIN, + END_MAX, + }; + IntervalPointType type; + }; + + std::vector interval_points_sorted_by_x_; + std::vector interval_points_sorted_by_y_; + + // Those two vectors are not strictly needed, we could instead iterate + // directly on the two vectors above, but the code would be much uglier. + struct PointsForCoordinate { + IntegerValue coordinate; + absl::Span points; + }; + std::vector grouped_intervals_sorted_by_x_; + std::vector grouped_intervals_sorted_by_y_; + + const std::vector& intervals_; + + IntegerValue minimum_energy_; + IntegerValue probe_area_; + int top_index_, bottom_index_, left_index_, right_index_; + absl::flat_hash_set ranges_touching_boundary_[4]; +}; + +// Monte-Carlo inspired heuristic to find a rectangles with an energy conflict: +// - start with a rectangle equals to the full bounding box of the elements; +// - shrink the rectangle by an edge to the next "interesting" value. Choose +// the edge randomly, but biased towards the change that increases the ratio +// area_inside / area_rectangle; +// - collect a result at every conflict; +// - stop when the rectangle is empty. +std::vector FindRectanglesWithEnergyConflictMC( + const std::vector& intervals, absl::BitGenRef random, + double temperature); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index fbded8c0e3..10f74b3406 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -2047,7 +2047,7 @@ void IntegerTrail::MergeReasonIntoInternal(std::vector* output) const { } // TODO(user): If this is called many time on the same variables, it could be -// made faster by using some caching mecanism. +// made faster by using some caching mechanism. absl::Span IntegerTrail::Reason(const Trail& trail, int trail_index) const { const int index = boolean_trail_index_to_integer_one_[trail_index]; @@ -2066,14 +2066,19 @@ absl::Span IntegerTrail::Reason(const Trail& trail, return *reason; } +void IntegerTrail::AppendNewBounds(std::vector* output) const { + return AppendNewBoundsFrom(vars_.size(), output); +} + // TODO(user): Implement a dense version if there is more trail entries // than variables! -void IntegerTrail::AppendNewBounds(std::vector* output) const { +void IntegerTrail::AppendNewBoundsFrom( + int base_index, std::vector* output) const { tmp_marked_.ClearAndResize(IntegerVariable(vars_.size())); // In order to push the best bound for each variable, we loop backward. - const int end = vars_.size(); - for (int i = integer_trail_.size(); --i >= end;) { + CHECK_GE(base_index, vars_.size()); + for (int i = integer_trail_.size(); --i >= base_index;) { const TrailEntry& entry = integer_trail_[i]; if (entry.var == kNoIntegerVariable) continue; if (tmp_marked_[entry.var]) continue; @@ -2296,7 +2301,8 @@ bool GenericLiteralWatcher::Propagate(Trail* trail) { void GenericLiteralWatcher::Untrail(const Trail& trail, int trail_index) { if (propagation_trail_index_ <= trail_index) { // Nothing to do since we found a conflict before Propagate() was called. - CHECK_EQ(propagation_trail_index_, trail_index); + CHECK_EQ(propagation_trail_index_, trail_index) + << " level " << trail.CurrentDecisionLevel(); return; } diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index b177267d00..e3a5a80f19 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -33,13 +33,9 @@ #include "absl/container/inlined_vector.h" #include "absl/log/check.h" #include "absl/strings/str_cat.h" -#include "absl/strings/string_view.h" #include "absl/types/span.h" -#include "ortools/base/hash.h" #include "ortools/base/logging.h" #include "ortools/base/strong_vector.h" -#include "ortools/base/types.h" -#include "ortools/graph/iterators.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" #include "ortools/sat/sat_parameters.pb.h" @@ -161,7 +157,7 @@ inline bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue* result) { // // Each time we create an IntegerVariable we also create its negation. This is // done like that so internally we only stores and deal with lower bound. The -// upper bound beeing the lower bound of the negated variable. +// upper bound being the lower bound of the negated variable. DEFINE_STRONG_INDEX_TYPE(IntegerVariable); const IntegerVariable kNoIntegerVariable(-1); inline IntegerVariable NegationOf(IntegerVariable i) { @@ -366,7 +362,7 @@ struct DebugSolution { std::vector proto_values; // This is filled from proto_values at load-time, and using the - // cp_model_mapping, we cache the solution of the integer-variabls that are + // cp_model_mapping, we cache the solution of the integer variables that are // mapped. Note that it is possible that not all integer variable are mapped. // // TODO(user): When this happen we should be able to infer the value of these @@ -586,7 +582,7 @@ class IntegerEncoder { for (const IntegerLiteral l : GetIntegerLiterals(lit)) { temp_associated_vars_.push_back(l.var); } - for (const auto [var, value] : GetEqualityLiterals(lit)) { + for (const auto& [var, value] : GetEqualityLiterals(lit)) { temp_associated_vars_.push_back(var); } return temp_associated_vars_; @@ -939,8 +935,8 @@ class IntegerTrail : public SatPropagator { // restrictive bound than the current one will have no effect. // // The reason for this "assignment" must be provided as: - // - A set of Literal currently beeing all false. - // - A set of IntegerLiteral currently beeing all true. + // - A set of Literal currently being all false. + // - A set of IntegerLiteral currently being all true. // // IMPORTANT: Notice the inversed sign in the literal reason. This is a bit // confusing but internally SAT use this direction for efficiency. @@ -1088,6 +1084,12 @@ class IntegerTrail : public SatPropagator { // propagations on the trail. void AppendNewBounds(std::vector* output) const; + // Inspects the trail and output all the non-level zero bounds from the base + // index (one per variables) to the output. The algo is sparse if there is + // only a few propagations on the trail. + void AppendNewBoundsFrom(int base_index, + std::vector* output) const; + // Returns the trail index < threshold of a TrailEntry about var. Returns -1 // if there is no such entry (at a positive decision level). This is basically // the trail index of the lower bound of var at the time. @@ -1103,8 +1105,8 @@ class IntegerTrail : public SatPropagator { IntegerVariable NextVariableToBranchOnInPropagationLoop() const; // If we had an incomplete propagation, it is important to fix all the - // variables and not relly on the propagation to do so. This is related to the - // InPropagationLoop() code above. + // variables and not really on the propagation to do so. This is related to + // the InPropagationLoop() code above. bool CurrentBranchHadAnIncompletePropagation(); IntegerVariable FirstUnassignedVariable() const; @@ -1209,7 +1211,7 @@ class IntegerTrail : public SatPropagator { // This is used by FindLowestTrailIndexThatExplainBound() and // FindTrailIndexOfVarBefore() to speed up the lookup. It keeps a trail index // for each variable that may or may not point to a TrailEntry regarding this - // variable. The validity of the index is verified before beeing used. + // variable. The validity of the index is verified before being used. // // The cache will only be updated with trail_index >= threshold. mutable int var_trail_index_cache_threshold_ = 0; @@ -1346,7 +1348,7 @@ class PropagatorInterface { // - At level zero, it will not contain any indices associated with literals // that were already fixed when the propagator was registered. Only the // indices of the literals modified after the registration will be present. - virtual bool IncrementalPropagate(const std::vector& watch_indices) { + virtual bool IncrementalPropagate(const std::vector& /*watch_indices*/) { LOG(FATAL) << "Not implemented."; return false; // Remove warning in Windows } diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index c0018b82a7..5bffd0c402 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -17,7 +17,6 @@ #include #include #include -#include #include #include #include @@ -27,7 +26,9 @@ #include "absl/meta/type_traits.h" #include "absl/random/distributions.h" #include "absl/strings/str_cat.h" +#include "absl/types/span.h" #include "ortools/base/logging.h" +#include "ortools/sat/clause.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/implied_bounds.h" @@ -741,7 +742,7 @@ std::function CumulativePrecedenceSearchHeuristic( added_demand[t] = demand_min; current_height += demand_min; } else if (first_skipped_task == -1) { - // We should have everyting needed here to add a new precedence. + // We should have everything needed here to add a new precedence. first_skipped_task = t; found = true; break; @@ -1249,10 +1250,15 @@ IntegerSearchHelper::IntegerSearchHelper(Model* model) prober_(model->GetOrCreate()), product_detector_(model->GetOrCreate()), time_limit_(model->GetOrCreate()), - pseudo_costs_(model->GetOrCreate()) { + pseudo_costs_(model->GetOrCreate()), + inprocessing_(model->GetOrCreate()) { // This is needed for recording the pseudo-costs. const ObjectiveDefinition* objective = model->Get(); if (objective != nullptr) objective_var_ = objective->objective_var; + + // This is needed for the first MaybeMinimizeByPropagation() not to trigger + // right away. + sat_solver_->ResetMinimizationByPropagationThreshold(); } bool IntegerSearchHelper::BeforeTakingDecision() { @@ -1267,21 +1273,27 @@ bool IntegerSearchHelper::BeforeTakingDecision() { } } - if (sat_solver_->CurrentDecisionLevel() == 0) { - auto* level_zero_callbacks = model_->GetOrCreate(); - for (const auto& cb : level_zero_callbacks->callbacks) { - if (!cb()) { - sat_solver_->NotifyThatModelIsUnsat(); - return false; - } - } + // The rest only trigger at level zero. + if (sat_solver_->CurrentDecisionLevel() != 0) return true; - if (parameters_.use_sat_inprocessing() && - !model_->GetOrCreate()->InprocessingRound()) { + auto* level_zero_callbacks = model_->GetOrCreate(); + for (const auto& cb : level_zero_callbacks->callbacks) { + if (!cb()) { sat_solver_->NotifyThatModelIsUnsat(); return false; } } + + if (parameters_.use_sat_inprocessing() && + !inprocessing_->InprocessingRound()) { + sat_solver_->NotifyThatModelIsUnsat(); + return false; + } + + if (!sat_solver_->MaybeMinimizeByPropagation()) { + return sat_solver_->UnsatStatus(); + } + return true; } @@ -1402,22 +1414,10 @@ SatSolver::Status IntegerSearchHelper::SolveIntegerProblem() { // how to fix. if (!sat_solver_->FinishPropagation()) return sat_solver_->UnsatStatus(); - // Used to trigger clause minimization via propagation. - // - // TODO(user): integrate this with BeforeTakingDecision() so that shared tree - // search worker can also use that. - int64_t num_restarts = 0; - int64_t next_minimization_num_restart = - parameters_.minimize_with_propagation_restart_period(); - if (next_minimization_num_restart < 0) { - next_minimization_num_restart = std::numeric_limits::max(); - } - // Main search loop. const int64_t old_num_conflicts = sat_solver_->num_failures(); const int64_t conflict_limit = parameters_.max_number_of_conflicts(); int64_t num_decisions_since_last_lp_record_ = 0; - int64_t num_decisions_without_probing = 0; while (!time_limit_->LimitReached() && (sat_solver_->num_failures() - old_num_conflicts < conflict_limit)) { // If needed, restart and switch decision_policy. @@ -1426,25 +1426,10 @@ SatSolver::Status IntegerSearchHelper::SolveIntegerProblem() { return sat_solver_->UnsatStatus(); } heuristics.policy_index = (heuristics.policy_index + 1) % num_policies; - ++num_restarts; } if (!BeforeTakingDecision()) return sat_solver_->UnsatStatus(); - // Clause minimization using propagation. - if (sat_solver_->CurrentDecisionLevel() == 0 && - num_restarts >= next_minimization_num_restart) { - next_minimization_num_restart = - num_restarts + parameters_.minimize_with_propagation_restart_period(); - sat_solver_->MinimizeSomeClauses( - parameters_.minimize_with_propagation_num_decisions()); - - // Note(user): In some corner cases, the function above might find a - // feasible assignment. I think it is okay to ignore this special case - // that should only happen on trivial problems and just reset the solver. - if (!sat_solver_->ResetToLevelZero()) return sat_solver_->UnsatStatus(); - } - LiteralIndex decision = kNoLiteralIndex; while (true) { if (sat_solver_->ModelIsUnsat()) return sat_solver_->UnsatStatus(); @@ -1464,29 +1449,6 @@ SatSolver::Status IntegerSearchHelper::SolveIntegerProblem() { continue; } } - - // Probing? - // - // TODO(user): Be smarter about what variables we probe, we can - // also do more than one. - if (decision != kNoLiteralIndex && - sat_solver_->CurrentDecisionLevel() == 0 && - parameters_.probing_period_at_root() > 0 && - ++num_decisions_without_probing >= - parameters_.probing_period_at_root()) { - num_decisions_without_probing = 0; - if (!prober_->ProbeOneVariable(Literal(decision).Variable())) { - return SatSolver::INFEASIBLE; - } - DCHECK_EQ(sat_solver_->CurrentDecisionLevel(), 0); - - // We need to check after the probing that the literal is not fixed, - // otherwise we just go to the next decision. - if (sat_solver_->Assignment().LiteralIsAssigned(Literal(decision))) { - decision = kNoLiteralIndex; - continue; - } - } break; } @@ -1586,22 +1548,39 @@ SatSolver::Status SolveIntegerProblemWithLazyEncoding(Model* model) { return ResetAndSolveIntegerProblem(/*assumptions=*/{}, model); } +#define RETURN_IF_NOT_FEASIBLE(test) \ + const SatSolver::Status status = (test); \ + if (status != SatSolver::FEASIBLE) return status; + ContinuousProber::ContinuousProber(const CpModelProto& model_proto, Model* model) : model_(model), sat_solver_(model->GetOrCreate()), time_limit_(model->GetOrCreate()), + binary_implication_graph_(model->GetOrCreate()), + literal_watchers_(model->GetOrCreate()), trail_(model->GetOrCreate()), integer_trail_(model->GetOrCreate()), encoder_(model->GetOrCreate()), + inprocessing_(model->GetOrCreate()), parameters_(*(model->GetOrCreate())), level_zero_callbacks_(model->GetOrCreate()), prober_(model->GetOrCreate()), shared_response_manager_(model->Mutable()), shared_bounds_manager_(model->Mutable()), + random_(model->GetOrCreate()), active_limit_(parameters_.shaving_search_deterministic_time()) { auto* mapping = model_->GetOrCreate(); absl::flat_hash_set visited; + + trail_index_at_start_of_iteration_ = trail_->Index(); + integer_trail_index_at_start_of_iteration_ = integer_trail_->Index(); + + // Build variable lists. + // TODO(user): Ideally, we should scan the internal model. But there can be + // a large blowup of variables during loading, which slows down the probing + // part. Using model variables is a good heuristic to select 'impactful' + // Boolean variables. for (int v = 0; v < model_proto.variables_size(); ++v) { if (mapping->IsBoolean(v)) { const BooleanVariable bool_var = mapping->Literal(v).Variable(); @@ -1615,10 +1594,10 @@ ContinuousProber::ContinuousProber(const CpModelProto& model_proto, int_vars_.push_back(var); } } + VLOG(2) << "Start continuous probing with " << bool_vars_.size() - << " Boolean variables, and " << int_vars_.size() - << " integer variables" - << ", deterministic time limit = " + << " Boolean variables, " << int_vars_.size() + << " integer variables, deterministic time limit = " << time_limit_->GetDeterministicLimit() << " on " << model_->Name(); } @@ -1634,15 +1613,12 @@ SatSolver::Status ContinuousProber::Probe() { if (!sat_solver_->ResetToLevelZero()) return SatSolver::INFEASIBLE; while (!time_limit_->LimitReached()) { - if (parameters_.minimize_with_propagation_restart_period() >= 0) { - sat_solver_->MinimizeSomeClauses( - parameters_.minimize_with_propagation_num_decisions()); + if (parameters_.use_sat_inprocessing() && + !inprocessing_->InprocessingRound()) { + sat_solver_->NotifyThatModelIsUnsat(); + return sat_solver_->UnsatStatus(); } - // Probe each Boolean variable at most once per loop. - probed_bool_vars_.clear(); - probed_literals_.clear(); - // Store current statistics to detect an iteration without any improvement. const int64_t initial_num_literals_fixed = prober_->num_new_literals_fixed(); @@ -1657,51 +1633,39 @@ SatSolver::Status ContinuousProber::Probe() { continue; } - if (!ImportFromSharedClasses()) { - return SatSolver::INFEASIBLE; - } - - if (time_limit_->LimitReached()) { - return SatSolver::LIMIT_REACHED; - } - - const BooleanVariable shave_lb = - encoder_ - ->GetOrCreateAssociatedLiteral(IntegerLiteral::LowerOrEqual( - int_var, integer_trail_->LowerBound(int_var))) - .Variable(); - const auto [_lb, lb_inserted] = probed_bool_vars_.insert(shave_lb); + const Literal shave_lb_literal = + encoder_->GetOrCreateAssociatedLiteral(IntegerLiteral::LowerOrEqual( + int_var, integer_trail_->LowerBound(int_var))); + const BooleanVariable shave_lb_var = shave_lb_literal.Variable(); + const auto [_lb, lb_inserted] = probed_bool_vars_.insert(shave_lb_var); if (lb_inserted) { - if (!prober_->ProbeOneVariable(shave_lb)) { + if (!prober_->ProbeOneVariable(shave_lb_var)) { return SatSolver::INFEASIBLE; } num_literals_probed_++; } - const BooleanVariable shave_ub = - encoder_ - ->GetOrCreateAssociatedLiteral(IntegerLiteral::GreaterOrEqual( - int_var, integer_trail_->UpperBound(int_var))) - .Variable(); - const auto [_ub, ub_inserted] = probed_bool_vars_.insert(shave_ub); + const Literal shave_ub_literal = + encoder_->GetOrCreateAssociatedLiteral(IntegerLiteral::GreaterOrEqual( + int_var, integer_trail_->UpperBound(int_var))); + const BooleanVariable shave_ub_var = shave_ub_literal.Variable(); + const auto [_ub, ub_inserted] = probed_bool_vars_.insert(shave_ub_var); if (ub_inserted) { - if (!prober_->ProbeOneVariable(shave_ub)) { + if (!prober_->ProbeOneVariable(shave_ub_var)) { return SatSolver::INFEASIBLE; } num_literals_probed_++; } - if (parameters_.use_shaving_in_probing_search()) { - const SatSolver::Status lb_status = - ShaveLiteral(Literal(shave_lb, true)); + if (use_shaving_) { + const SatSolver::Status lb_status = ShaveLiteral(shave_lb_literal); if (ReportStatus(lb_status)) return lb_status; - const SatSolver::Status ub_status = - ShaveLiteral(Literal(shave_ub, true)); + const SatSolver::Status ub_status = ShaveLiteral(shave_ub_literal); if (ReportStatus(ub_status)) return ub_status; } - LogStatistics(); + RETURN_IF_NOT_FEASIBLE(PeriodicSyncAndCheck()); } // Probe Boolean variables from the model. @@ -1710,24 +1674,16 @@ SatSolver::Status ContinuousProber::Probe() { if (sat_solver_->Assignment().VariableIsAssigned(bool_var)) continue; - if (!ImportFromSharedClasses()) { + const auto [_, inserted] = probed_bool_vars_.insert(bool_var); + if (!inserted) continue; + + if (!prober_->ProbeOneVariable(bool_var)) { return SatSolver::INFEASIBLE; } - - if (time_limit_->LimitReached()) { - return SatSolver::LIMIT_REACHED; - } - - const auto [_, inserted] = probed_bool_vars_.insert(bool_var); - if (inserted) { - if (!prober_->ProbeOneVariable(bool_var)) { - return SatSolver::INFEASIBLE; - } - num_literals_probed_++; - } + num_literals_probed_++; const Literal literal(bool_var, true); - if (parameters_.use_shaving_in_probing_search() && + if (use_shaving_ && !sat_solver_->Assignment().LiteralIsAssigned(literal)) { const SatSolver::Status true_status = ShaveLiteral(literal); if (ReportStatus(true_status)) return true_status; @@ -1737,39 +1693,198 @@ SatSolver::Status ContinuousProber::Probe() { if (ReportStatus(false_status)) return false_status; } - LogStatistics(); + RETURN_IF_NOT_FEASIBLE(PeriodicSyncAndCheck()); + } + + if (parameters_.use_extended_probing()) { + const auto at_least_one_literal_is_true = + [this](absl::Span literals) { + for (const Literal literal : literals) { + if (sat_solver_->Assignment().LiteralIsTrue(literal)) { + return true; + } + } + return false; + }; + + // Probe clauses of the SAT model. + for (;;) { + const SatClause* clause = literal_watchers_->NextClauseToProbe(); + if (clause == nullptr) break; + if (at_least_one_literal_is_true(clause->AsSpan())) continue; + + tmp_dnf_.clear(); + for (const Literal literal : clause->AsSpan()) { + if (sat_solver_->Assignment().LiteralIsAssigned(literal)) continue; + tmp_dnf_.push_back({literal}); + } + ++num_at_least_one_probed_; + if (!prober_->ProbeDnf("at_least_one", tmp_dnf_)) { + return SatSolver::INFEASIBLE; + } + + RETURN_IF_NOT_FEASIBLE(PeriodicSyncAndCheck()); + } + + // Probe at_most_ones of the SAT model. + for (;;) { + const absl::Span at_most_one = + binary_implication_graph_->NextAtMostOne(); + if (at_most_one.empty()) break; + if (at_least_one_literal_is_true(at_most_one)) continue; + + tmp_dnf_.clear(); + tmp_literals_.clear(); + for (const Literal literal : at_most_one) { + if (sat_solver_->Assignment().LiteralIsAssigned(literal)) continue; + tmp_dnf_.push_back({literal}); + tmp_literals_.push_back(literal.Negated()); + } + tmp_dnf_.push_back(tmp_literals_); + ++num_at_most_one_probed_; + if (!prober_->ProbeDnf("at_most_one", tmp_dnf_)) { + return SatSolver::INFEASIBLE; + } + + RETURN_IF_NOT_FEASIBLE(PeriodicSyncAndCheck()); + } + + // Probe combinations of Booleans variables. + // TODO(user): Probe until something is imported, or learned. + if (bool_vars_.size() < 1000) { + for (; current_bv1_ + 1 < bool_vars_.size(); ++current_bv1_) { + const BooleanVariable bv1 = bool_vars_[current_bv1_]; + if (sat_solver_->Assignment().VariableIsAssigned(bv1)) continue; + current_bv2_ = std::max(current_bv1_ + 1, current_bv2_); + for (; current_bv2_ < bool_vars_.size(); ++current_bv2_) { + const BooleanVariable& bv2 = bool_vars_[current_bv2_]; + if (sat_solver_->Assignment().VariableIsAssigned(bv2)) continue; + if (!prober_->ProbeDnf( + "pair_of_bool_vars", + {{Literal(bv1, true), Literal(bv2, true)}, + {Literal(bv1, true), Literal(bv2, false)}, + {Literal(bv1, false), Literal(bv2, true)}, + {Literal(bv1, false), Literal(bv2, false)}})) { + return SatSolver::INFEASIBLE; + } + RETURN_IF_NOT_FEASIBLE(PeriodicSyncAndCheck()); + } + current_bv2_ = 0; + } + } else { + for (; random_pair_of_bool_vars_probed_ < 5000; + ++random_pair_of_bool_vars_probed_) { + const BooleanVariable bv1 = + bool_vars_[absl::Uniform(*random_, 0, bool_vars_.size())]; + if (sat_solver_->Assignment().VariableIsAssigned(bv1)) continue; + const BooleanVariable bv2 = + bool_vars_[absl::Uniform(*random_, 0, bool_vars_.size())]; + if (sat_solver_->Assignment().VariableIsAssigned(bv2) || bv1 == bv2) { + continue; + } + if (!prober_->ProbeDnf( + "rnd_pair_of_bool_vars", + {{Literal(bv1, true), Literal(bv2, true)}, + {Literal(bv1, true), Literal(bv2, false)}, + {Literal(bv1, false), Literal(bv2, true)}, + {Literal(bv1, false), Literal(bv2, false)}})) { + return SatSolver::INFEASIBLE; + } + + RETURN_IF_NOT_FEASIBLE(PeriodicSyncAndCheck()); + } + } } // Adjust the active_limit. { - const double deterministic_time = - parameters_.shaving_search_deterministic_time(); - const bool something_has_been_detected = - num_bounds_shaved_ != initial_num_bounds_shaved || - prober_->num_new_literals_fixed() != initial_num_literals_fixed; - if (something_has_been_detected) { // Reset the limit. - active_limit_ = deterministic_time; - } else if (active_limit_ < 25 * deterministic_time) { // Bump the limit. - active_limit_ += deterministic_time; + if (use_shaving_) { + const double deterministic_time = + parameters_.shaving_search_deterministic_time(); + const bool something_has_been_detected = + num_bounds_shaved_ != initial_num_bounds_shaved || + prober_->num_new_literals_fixed() != initial_num_literals_fixed; + if (something_has_been_detected) { // Reset the limit. + active_limit_ = deterministic_time; + } else if (active_limit_ < + 25 * deterministic_time) { // Bump the limit. + active_limit_ += deterministic_time; + } } } + // Reset all counters. ++iteration_; current_bool_var_ = 0; current_int_var_ = 0; + current_bv1_ = 0; + current_bv2_ = 1; + random_pair_of_bool_vars_probed_ = 0; + binary_implication_graph_->ResetAtMostOneIterator(); + literal_watchers_->ResetToProbeIndex(); + probed_bool_vars_.clear(); + probed_literals_.clear(); + + const int new_trail_index = trail_->Index(); + const int new_integer_trail_index = integer_trail_->Index(); + + // Update the use_shaving_ parameter. + // TODO(user): Currently, the heuristics is that we alternate shaving and + // not shaving, unless use_shaving_in_probing_search is false. + use_shaving_ = + parameters_.use_shaving_in_probing_search() ? !use_shaving_ : false; + trail_index_at_start_of_iteration_ = new_trail_index; + integer_trail_index_at_start_of_iteration_ = new_integer_trail_index; + + // Remove fixed Boolean variables. + int size = 0; + for (int i = 0; i < bool_vars_.size(); ++i) { + if (!sat_solver_->Assignment().VariableIsAssigned(bool_vars_[i])) { + bool_vars_[size++] = bool_vars_[i]; + } + } + bool_vars_.resize(size); + + // Remove fixed integer variables. + size = 0; + for (int i = 0; i < int_vars_.size(); ++i) { + if (!integer_trail_->IsFixed(int_vars_[i])) { + int_vars_[size++] = int_vars_[i]; + } + } + int_vars_.resize(size); } return SatSolver::LIMIT_REACHED; } -bool ContinuousProber::ImportFromSharedClasses() { - if (!sat_solver_->ResetToLevelZero()) return false; - for (const auto& cb : level_zero_callbacks_->callbacks) { - if (!cb()) { - sat_solver_->NotifyThatModelIsUnsat(); - return false; +#undef RETURN_IF_NOT_FEASIBLE + +SatSolver::Status ContinuousProber::PeriodicSyncAndCheck() { + // Check limit. + if (--num_probes_remaining_ <= 0) { + num_probes_remaining_ = kTestLimitPeriod; + if (time_limit_->LimitReached()) return SatSolver::LIMIT_REACHED; + } + + // Log the state of the prober. + if (--num_logs_remaining_ <= 0) { + num_logs_remaining_ = kLogPeriod; + LogStatistics(); + } + + // Sync with shared storage. + if (--num_syncs_remaining_ <= 0) { + num_syncs_remaining_ = kSyncPeriod; + if (!sat_solver_->ResetToLevelZero()) return SatSolver::INFEASIBLE; + for (const auto& cb : level_zero_callbacks_->callbacks) { + if (!cb()) { + sat_solver_->NotifyThatModelIsUnsat(); + return SatSolver::INFEASIBLE; + } } } - return true; + + return SatSolver::FEASIBLE; } SatSolver::Status ContinuousProber::ShaveLiteral(Literal literal) { @@ -1812,11 +1927,16 @@ void ContinuousProber::LogStatistics() { "Probe", absl::StrCat( " (iterations=", iteration_, + " linearization_level=", parameters_.linearization_level(), + " shaving=", use_shaving_, " active_bool_vars=", bool_vars_.size(), + " active_int_vars=", integer_trail_->NumIntegerVariables(), " literals fixed/probed=", prober_->num_new_literals_fixed(), "/", num_literals_probed_, " bounds shaved/tried=", num_bounds_shaved_, "/", num_bounds_tried_, " new_integer_bounds=", shared_bounds_manager_->NumBoundsExported("probing"), - ", new_binary_clause=", prober_->num_new_binary_clauses(), ")")); + " new_binary_clause=", prober_->num_new_binary_clauses(), + " num_at_least_one_probed=", num_at_least_one_probed_, + " num_at_most_one_probed=", num_at_most_one_probed_, ")")); } } diff --git a/ortools/sat/integer_search.h b/ortools/sat/integer_search.h index 5e66e801b8..1017480865 100644 --- a/ortools/sat/integer_search.h +++ b/ortools/sat/integer_search.h @@ -29,7 +29,7 @@ #include #include "absl/container/flat_hash_set.h" -#include "absl/time/time.h" +#include "ortools/sat/clause.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/implied_bounds.h" #include "ortools/sat/integer.h" @@ -37,9 +37,11 @@ #include "ortools/sat/probing.h" #include "ortools/sat/pseudo_costs.h" #include "ortools/sat/sat_base.h" +#include "ortools/sat/sat_inprocessing.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/sat_solver.h" #include "ortools/sat/synchronization.h" +#include "ortools/sat/util.h" #include "ortools/util/strong_integers.h" #include "ortools/util/time_limit.h" @@ -254,12 +256,6 @@ std::vector> CompleteHeuristics( incomplete_heuristics, const std::function& completion_heuristic); -// Specialized search that will continuously probe Boolean variables and bounds -// of integer variables. -SatSolver::Status ContinuousProbing( - const std::vector& bool_vars, - const std::vector& int_vars, Model* model); - // An helper class to share the code used by the different kind of search. class IntegerSearchHelper { public: @@ -309,6 +305,7 @@ class IntegerSearchHelper { ProductDetector* product_detector_; TimeLimit* time_limit_; PseudoCosts* pseudo_costs_; + Inprocessing* inprocessing_; IntegerVariable objective_var_ = kNoIntegerVariable; bool must_process_conflict_ = false; @@ -331,10 +328,14 @@ class ContinuousProber { SatSolver::Status Probe(); private: - bool ImportFromSharedClasses(); + static const int kTestLimitPeriod = 50; + static const int kLogPeriod = 1000; + static const int kSyncPeriod = 50; + SatSolver::Status ShaveLiteral(Literal literal); bool ReportStatus(SatSolver::Status status); void LogStatistics(); + SatSolver::Status PeriodicSyncAndCheck(); // Variables to probe. std::vector bool_vars_; @@ -344,19 +345,35 @@ class ContinuousProber { Model* model_; SatSolver* sat_solver_; TimeLimit* time_limit_; + BinaryImplicationGraph* binary_implication_graph_; + LiteralWatchers* literal_watchers_; Trail* trail_; IntegerTrail* integer_trail_; IntegerEncoder* encoder_; + Inprocessing* inprocessing_; const SatParameters parameters_; LevelZeroCallbackHelper* level_zero_callbacks_; Prober* prober_; SharedResponseManager* shared_response_manager_; SharedBoundsManager* shared_bounds_manager_; + ModelRandomGenerator* random_; // Statistics. int64_t num_literals_probed_ = 0; int64_t num_bounds_shaved_ = 0; int64_t num_bounds_tried_ = 0; + int64_t num_at_least_one_probed_ = 0; + int64_t num_at_most_one_probed_ = 0; + + // Period counters; + int num_logs_remaining_ = 0; + int num_syncs_remaining_ = 0; + int num_probes_remaining_ = 0; + + // Shaving management. + bool use_shaving_ = false; + int trail_index_at_start_of_iteration_ = 0; + int integer_trail_index_at_start_of_iteration_ = 0; // Current state of the probe. double active_limit_; @@ -366,6 +383,11 @@ class ContinuousProber { int iteration_ = 1; int current_int_var_ = 0; int current_bool_var_ = 0; + int current_bv1_ = 0; + int current_bv2_ = 0; + int random_pair_of_bool_vars_probed_ = 0; + std::vector> tmp_dnf_; + std::vector tmp_literals_; }; } // namespace sat diff --git a/ortools/sat/intervals.h b/ortools/sat/intervals.h index 21dfeeccb4..ebb3850a0d 100644 --- a/ortools/sat/intervals.h +++ b/ortools/sat/intervals.h @@ -333,6 +333,13 @@ class SchedulingConstraintHelper : public PropagatorInterface, IntegerValue StartMax(int t) const { return -cached_negated_start_max_[t]; } IntegerValue EndMax(int t) const { return -cached_negated_end_max_[t]; } + IntegerValue LevelZeroStartMin(int t) const { + return integer_trail_->LevelZeroLowerBound(starts_[t]); + } + IntegerValue LevelZeroStartMax(int t) const { + return integer_trail_->LevelZeroUpperBound(starts_[t]); + } + // In the presence of tasks with a variable size, we do not necessarily // have start_min + size_min = end_min, we can instead have a situation // like: diff --git a/ortools/sat/parameters_validation.cc b/ortools/sat/parameters_validation.cc index c6761c3a30..19646910e7 100644 --- a/ortools/sat/parameters_validation.cc +++ b/ortools/sat/parameters_validation.cc @@ -82,6 +82,7 @@ std::string ValidateParameters(const SatParameters& params) { TEST_IS_FINITE(cut_max_active_count_value); TEST_IS_FINITE(cut_active_count_decay); TEST_IS_FINITE(shaving_search_deterministic_time); + TEST_IS_FINITE(minimize_with_propagation_ratio); TEST_IS_FINITE(mip_max_bound); TEST_IS_FINITE(mip_wanted_precision); TEST_IS_FINITE(mip_check_precision); @@ -89,6 +90,8 @@ std::string ValidateParameters(const SatParameters& params) { TEST_IS_FINITE(mip_drop_tolerance); TEST_IS_FINITE(shared_tree_worker_objective_split_probability); + TEST_POSITIVE(at_most_one_max_expansion_size); + TEST_NOT_NAN(max_time_in_seconds); TEST_NOT_NAN(max_deterministic_time); diff --git a/ortools/sat/pb_constraint.h b/ortools/sat/pb_constraint.h index 6afe72b28e..7eda79ddc1 100644 --- a/ortools/sat/pb_constraint.h +++ b/ortools/sat/pb_constraint.h @@ -27,7 +27,6 @@ #include "absl/strings/string_view.h" #include "absl/types/span.h" #include "ortools/base/logging.h" -#include "ortools/base/macros.h" #include "ortools/base/strong_vector.h" #include "ortools/base/types.h" #include "ortools/sat/model.h" @@ -153,6 +152,11 @@ class CanonicalBooleanLinearProblem { public: CanonicalBooleanLinearProblem() = default; + // This type is neither copyable nor movable. + CanonicalBooleanLinearProblem(const CanonicalBooleanLinearProblem&) = delete; + CanonicalBooleanLinearProblem& operator=( + const CanonicalBooleanLinearProblem&) = delete; + // Adds a new constraint to the problem. The bounds are inclusive. // Returns false in case of a possible overflow or if the constraint is // never satisfiable. @@ -175,7 +179,6 @@ class CanonicalBooleanLinearProblem { std::vector rhs_; std::vector> constraints_; - DISALLOW_COPY_AND_ASSIGN(CanonicalBooleanLinearProblem); }; // Encode a constraint sum term <= rhs, where each term is a positive @@ -537,6 +540,10 @@ class PbConstraints : public SatPropagator { num_threshold_updates_(0) { model->GetOrCreate()->RegisterPropagator(this); } + + // This type is neither copyable nor movable. + PbConstraints(const PbConstraints&) = delete; + PbConstraints& operator=(const PbConstraints&) = delete; ~PbConstraints() override { IF_STATS_ENABLED({ LOG(INFO) << stats_.StatString(); @@ -686,7 +693,6 @@ class PbConstraints : public SatPropagator { int64_t num_constraint_lookups_; int64_t num_inspected_constraint_literals_; int64_t num_threshold_updates_; - DISALLOW_COPY_AND_ASSIGN(PbConstraints); }; // Boolean linear constraints can propagate a lot of literals at the same time. @@ -701,6 +707,12 @@ class VariableWithSameReasonIdentifier { explicit VariableWithSameReasonIdentifier(const Trail& trail) : trail_(trail) {} + // This type is neither copyable nor movable. + VariableWithSameReasonIdentifier(const VariableWithSameReasonIdentifier&) = + delete; + VariableWithSameReasonIdentifier& operator=( + const VariableWithSameReasonIdentifier&) = delete; + void Resize(int num_variables) { first_variable_.resize(num_variables); seen_.ClearAndResize(BooleanVariable(num_variables)); @@ -727,8 +739,6 @@ class VariableWithSameReasonIdentifier { const Trail& trail_; absl::StrongVector first_variable_; SparseBitset seen_; - - DISALLOW_COPY_AND_ASSIGN(VariableWithSameReasonIdentifier); }; } // namespace sat diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index 9d4f677501..73f46129fa 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -735,6 +735,8 @@ bool PresolveContext::PropagateAffineRelation(int ref) { bool PresolveContext::PropagateAffineRelation(int ref, int rep, int64_t coeff, int64_t offset) { + DCHECK(!DomainIsEmpty(ref)); + DCHECK(!DomainIsEmpty(rep)); if (!RefIsPositive(rep)) { rep = NegatedRef(rep); coeff = -coeff; @@ -1299,7 +1301,10 @@ void PresolveContext::InsertVarValueEncodingInternal(int literal, int var, bool PresolveContext::InsertHalfVarValueEncoding(int literal, int var, int64_t value, bool imply_eq) { if (is_unsat_) return false; - CHECK(RefIsPositive(var)); + DCHECK(RefIsPositive(var)); + if (!CanonicalizeEncoding(&var, &value) || !DomainOf(var).Contains(value)) { + return SetLiteralToFalse(literal); + } // Creates the linking sets on demand. // Insert the enforcement literal in the half encoding map. @@ -1338,7 +1343,7 @@ bool PresolveContext::CanonicalizeEncoding(int* ref, int64_t* value) { bool PresolveContext::InsertVarValueEncoding(int literal, int ref, int64_t value) { - if (!CanonicalizeEncoding(&ref, &value)) { + if (!CanonicalizeEncoding(&ref, &value) || !DomainOf(ref).Contains(value)) { return SetLiteralToFalse(literal); } literal = GetLiteralRepresentative(literal); @@ -1348,7 +1353,9 @@ bool PresolveContext::InsertVarValueEncoding(int literal, int ref, bool PresolveContext::StoreLiteralImpliesVarEqValue(int literal, int var, int64_t value) { - if (!CanonicalizeEncoding(&var, &value)) return false; + if (!CanonicalizeEncoding(&var, &value) || !DomainOf(var).Contains(value)) { + return SetLiteralToFalse(literal); + } literal = GetLiteralRepresentative(literal); return InsertHalfVarValueEncoding(literal, var, value, /*imply_eq=*/true); } diff --git a/ortools/sat/probing.cc b/ortools/sat/probing.cc index 95c682e787..382488ffb8 100644 --- a/ortools/sat/probing.cc +++ b/ortools/sat/probing.cc @@ -19,8 +19,11 @@ #include #include +#include "absl/container/btree_map.h" +#include "absl/container/btree_set.h" #include "absl/container/inlined_vector.h" #include "absl/log/check.h" +#include "absl/strings/string_view.h" #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/strong_vector.h" @@ -298,6 +301,122 @@ bool Prober::ProbeBooleanVariables( return true; } +bool Prober::ProbeDnf(absl::string_view name, + const std::vector>& dnf) { + if (dnf.size() <= 1) return true; + + // Reset the solver in case it was already used. + if (!sat_solver_->ResetToLevelZero()) return false; + + always_propagated_bounds_.clear(); + always_propagated_literals_.clear(); + int num_valid_conjunctions = 0; + for (const std::vector& conjunction : dnf) { + if (!sat_solver_->ResetToLevelZero()) return false; + if (num_valid_conjunctions > 0 && always_propagated_bounds_.empty() && + always_propagated_literals_.empty()) { + // We can exit safely as nothing will be propagated. + return true; + } + + bool conjunction_is_valid = true; + int num_literals_enqueued = 0; + const int root_trail_index = trail_.Index(); + const int root_integer_trail_index = integer_trail_->Index(); + for (const Literal& lit : conjunction) { + if (assignment_.LiteralIsAssigned(lit)) { + if (assignment_.LiteralIsTrue(lit)) continue; + conjunction_is_valid = false; + break; + } + ++num_literals_enqueued; + const int decision_level_before_enqueue = + sat_solver_->CurrentDecisionLevel(); + sat_solver_->EnqueueDecisionAndBackjumpOnConflict(lit); + sat_solver_->AdvanceDeterministicTime(time_limit_); + const int decision_level_after_enqueue = + sat_solver_->CurrentDecisionLevel(); + ++num_decisions_; + + if (sat_solver_->ModelIsUnsat()) return false; + // If the literal has been pushed without any conflict, the level should + // have been increased. + if (decision_level_after_enqueue <= decision_level_before_enqueue) { + conjunction_is_valid = false; + break; + } + // TODO(user): Can we use the callback_? + } + if (conjunction_is_valid) { + ++num_valid_conjunctions; + } else { + continue; + } + + // Process propagated literals. + new_propagated_literals_.clear(); + for (int i = root_trail_index; i < trail_.Index(); ++i) { + const LiteralIndex literal_index = trail_[i].Index(); + if (num_valid_conjunctions == 1 || + always_propagated_literals_.contains(literal_index)) { + new_propagated_literals_.insert(literal_index); + } + } + std::swap(new_propagated_literals_, always_propagated_literals_); + + // Process propagated integer bounds. + new_integer_bounds_.clear(); + integer_trail_->AppendNewBoundsFrom(root_integer_trail_index, + &new_integer_bounds_); + new_propagated_bounds_.clear(); + for (const IntegerLiteral entry : new_integer_bounds_) { + const auto it = always_propagated_bounds_.find(entry.var); + if (num_valid_conjunctions == 1) { // First loop. + new_propagated_bounds_[entry.var] = entry.bound; + } else if (it != always_propagated_bounds_.end()) { + new_propagated_bounds_[entry.var] = std::min(entry.bound, it->second); + } + } + std::swap(new_propagated_bounds_, always_propagated_bounds_); + } + + if (!sat_solver_->ResetToLevelZero()) return false; + // Fix literals implied by the dnf. + const int previous_num_literals_fixed = num_new_literals_fixed_; + for (const LiteralIndex literal_index : always_propagated_literals_) { + const Literal lit(literal_index); + if (assignment_.LiteralIsTrue(lit)) continue; + ++num_new_literals_fixed_; + if (!sat_solver_->AddUnitClause(lit)) return false; + } + + // Fix integer bounds implied by the dnf. + int previous_num_integer_bounds = num_new_integer_bounds_; + for (const auto& [var, bound] : always_propagated_bounds_) { + if (bound > integer_trail_->LowerBound(var)) { + ++num_new_integer_bounds_; + if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(var, bound), + {}, {})) { + return false; + } + } + } + + if (!sat_solver_->FinishPropagation()) return false; + + if (num_new_integer_bounds_ > previous_num_integer_bounds || + num_new_literals_fixed_ > previous_num_literals_fixed) { + VLOG(1) << "ProbeDnf(" << name << ", num_fixed_literals=" + << num_new_literals_fixed_ - previous_num_literals_fixed + << ", num_fixed_integer_bounds=" + << num_new_integer_bounds_ - previous_num_integer_bounds + << ", num_valid_conjunctions=" << num_valid_conjunctions << "/" + << dnf.size() << ")"; + } + + return true; +} + bool LookForTrivialSatSolution(double deterministic_time_limit, Model* model, SolverLogger* logger) { WallTimer wall_timer; @@ -646,7 +765,7 @@ bool FailedLiteralProbingRound(ProbingOptions options, Model* model) { if (subsumed) { ++num_new_subsumed; ++num_new_binary; - implication_graph->AddBinaryClause(last_decision.Negated(), l); + CHECK(implication_graph->AddBinaryClause(last_decision.Negated(), l)); const int trail_index = trail.Info(l.Variable()).trail_index; int test = 0; @@ -683,7 +802,7 @@ bool FailedLiteralProbingRound(ProbingOptions options, Model* model) { // this clause subsume the reason. if (!subsumed && trail.AssignmentType(l.Variable()) != id) { ++num_new_binary; - implication_graph->AddBinaryClause(last_decision.Negated(), l); + CHECK(implication_graph->AddBinaryClause(last_decision.Negated(), l)); } } else { // If we don't extract binary, we don't need to explore any of @@ -705,7 +824,7 @@ bool FailedLiteralProbingRound(ProbingOptions options, Model* model) { for (const auto& w : clause_manager->WatcherListOnFalse(last_decision.Negated())) { if (assignment.LiteralIsTrue(w.blocking_literal)) { - if (w.clause->empty()) continue; + if (w.clause->IsRemoved()) continue; CHECK_NE(w.blocking_literal, last_decision.Negated()); // Add the binary clause if needed. Note that we change the reason @@ -713,7 +832,7 @@ bool FailedLiteralProbingRound(ProbingOptions options, Model* model) { // // Tricky: while last_decision would be a valid reason, we need a // reason that was assigned before this literal, so we use the - // decision at the level where this literal was assigne which is an + // decision at the level where this literal was assigned which is an // even better reasony. Maybe it is just better to change all the // reason above to a binary one so we don't have an issue here. if (trail.AssignmentType(w.blocking_literal.Variable()) != id) { @@ -722,8 +841,8 @@ bool FailedLiteralProbingRound(ProbingOptions options, Model* model) { const auto& info = trail.Info(w.blocking_literal.Variable()); if (info.level > 0) { ++num_new_binary; - implication_graph->AddBinaryClause(last_decision.Negated(), - w.blocking_literal); + CHECK(implication_graph->AddBinaryClause(last_decision.Negated(), + w.blocking_literal)); const Literal d = sat_solver->Decisions()[info.level - 1].literal; if (d != w.blocking_literal) { diff --git a/ortools/sat/probing.h b/ortools/sat/probing.h index 349faa9012..59f71a443f 100644 --- a/ortools/sat/probing.h +++ b/ortools/sat/probing.h @@ -19,7 +19,10 @@ #include #include +#include "absl/container/btree_map.h" +#include "absl/container/btree_set.h" #include "absl/strings/str_cat.h" +#include "absl/strings/string_view.h" #include "absl/types/span.h" #include "ortools/sat/clause.h" #include "ortools/sat/implied_bounds.h" @@ -27,7 +30,6 @@ #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" #include "ortools/sat/sat_solver.h" -#include "ortools/sat/util.h" #include "ortools/util/bitset.h" #include "ortools/util/logging.h" #include "ortools/util/time_limit.h" @@ -81,6 +83,12 @@ class Prober { bool ProbeOneVariable(BooleanVariable b); + // Probes the given problem DNF (disjunction of conjunctions). Since one of + // the conjunction must be true, we might be able to fix literal or improve + // integer bounds if all conjunction propagate the same thing. + bool ProbeDnf(absl::string_view name, + const std::vector>& dnf); + // Statistics. // They are reset each time ProbleBooleanVariables() is called. // Note however that we do not reset them on a call to ProbeOneVariable(). @@ -117,6 +125,10 @@ class Prober { std::vector to_fix_at_true_; std::vector new_integer_bounds_; std::vector> new_binary_clauses_; + absl::btree_set new_propagated_literals_; + absl::btree_set always_propagated_literals_; + absl::btree_map new_propagated_bounds_; + absl::btree_map always_propagated_bounds_; // Probing statistics. int num_decisions_ = 0; @@ -199,7 +211,7 @@ struct ProbingOptions { // resolution have been performed may need to do more work. bool use_tree_look = true; - // There is two sligthly different implementation of the tree-look algo. + // There is two slightly different implementation of the tree-look algo. // // TODO(user): Decide which one is better, currently the difference seems // small but the queue seems slightly faster. @@ -208,7 +220,7 @@ struct ProbingOptions { // If we detect as we probe that a new binary clause subsumes one of the // non-binary clause, we will replace the long clause by the binary one. This // is orthogonal to the extract_binary_clauses parameters which will add all - // binary clauses but not neceassirly check for subsumption. + // binary clauses but not necessarily check for subsumption. bool subsume_with_binary_clause = true; // We assume this is also true if --v 1 is activated. diff --git a/ortools/sat/sat_base.h b/ortools/sat/sat_base.h index 5eedef7056..6072d2e37f 100644 --- a/ortools/sat/sat_base.h +++ b/ortools/sat/sat_base.h @@ -20,7 +20,6 @@ #include #include #include -#include #include #include #include @@ -33,7 +32,6 @@ #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/strong_vector.h" -#include "ortools/base/types.h" #include "ortools/sat/model.h" #include "ortools/util/bitset.h" #include "ortools/util/strong_integers.h" @@ -120,6 +118,11 @@ inline std::ostream& operator<<(std::ostream& os, Literal literal) { return os; } +template +void AbslStringify(Sink& sink, Literal arg) { + absl::Format(&sink, "%s", arg.DebugString()); +} + inline std::ostream& operator<<(std::ostream& os, absl::Span literals) { os << "["; @@ -414,7 +417,7 @@ class Trail { // Returns the last conflict. absl::Span FailingClause() const { if (DEBUG_MODE && debug_checker_ != nullptr) { - debug_checker_(conflict_); + CHECK(debug_checker_(conflict_)); } return conflict_; } @@ -443,7 +446,7 @@ class Trail { } // Print the current literals on the trail. - std::string DebugString() { + std::string DebugString() const { std::string result; for (int i = 0; i < current_info_.trail_index; ++i) { if (!result.empty()) result += " "; @@ -555,7 +558,7 @@ class SatPropagator { // can use trail_.GetEmptyVectorToStoreReason() if it doesn't have a memory // location that already contains the reason. virtual absl::Span Reason(const Trail& trail, - int trail_index) const { + int /*trail_index*/) const { LOG(FATAL) << "Not implemented."; return {}; } @@ -659,7 +662,7 @@ inline absl::Span Trail::Reason(BooleanVariable var) const { std::vector clause; clause.assign(reasons_[var].begin(), reasons_[var].end()); clause.push_back(assignment_.GetTrueLiteralForAssignedVariable(var)); - debug_checker_(clause); + CHECK(debug_checker_(clause)); } return reasons_[var]; } @@ -679,7 +682,7 @@ inline absl::Span Trail::Reason(BooleanVariable var) const { std::vector clause; clause.assign(reasons_[var].begin(), reasons_[var].end()); clause.push_back(assignment_.GetTrueLiteralForAssignedVariable(var)); - debug_checker_(clause); + CHECK(debug_checker_(clause)); } return reasons_[var]; } diff --git a/ortools/sat/sat_inprocessing.cc b/ortools/sat/sat_inprocessing.cc index 8b136dcdd9..fe3483c067 100644 --- a/ortools/sat/sat_inprocessing.cc +++ b/ortools/sat/sat_inprocessing.cc @@ -70,7 +70,7 @@ bool Inprocessing::PresolveLoop(SatPresolveOptions options) { const bool log_round_info = VLOG_IS_ON(1); // We currently do the transformations in a given order and restart each time - // we did something to make sure that the earlier step cannot srengthen more. + // we did something to make sure that the earlier step cannot strengthen more. // This might not be the best, but it is really good during development phase // to make sure each individual functions is as incremental and as fast as // possible. @@ -94,6 +94,14 @@ bool Inprocessing::PresolveLoop(SatPresolveOptions options) { // seems better to do just one loop instead of two over all clauses. Because // of memory access. it isn't that clear though. RETURN_IF_FALSE(RemoveFixedAndEquivalentVariables(log_round_info)); + + // IMPORTANT: Since we only run this on pure sat problem, we can just + // get rid of equivalent variable right away and do not need to keep them + // in the implication_graph_ for propagation. + // + // This is needed for the correctness of the bounded variable elimination. + implication_graph_->RemoveAllRedundantVariables(&postsolve_->clauses); + RETURN_IF_FALSE(stamping_simplifier_->DoOneRound(log_round_info)); // We wait for the fix-point to be reached before doing the other @@ -113,6 +121,9 @@ bool Inprocessing::PresolveLoop(SatPresolveOptions options) { // clause graph twice. It might make sense to reach the BCE fix point which // is unique before each variable elimination. blocked_clause_simplifier_->DoOneRound(log_round_info); + + // TODO(user): this break some binary graph invariant. Fix! + RETURN_IF_FALSE(RemoveFixedAndEquivalentVariables(log_round_info)); RETURN_IF_FALSE(bounded_variable_elimination_->DoOneRound(log_round_info)); RETURN_IF_FALSE(LevelZeroPropagate()); @@ -157,8 +168,8 @@ bool Inprocessing::InprocessingRound() { WallTimer wall_timer; wall_timer.Start(); - const bool log_info = true || VLOG_IS_ON(1); - const bool log_round_info = VLOG_IS_ON(1); + const bool log_info = VLOG_IS_ON(1); + const bool log_round_info = VLOG_IS_ON(2); // Mainly useful for development. double probing_time = 0.0; @@ -193,31 +204,46 @@ bool Inprocessing::InprocessingRound() { RETURN_IF_FALSE(stamping_simplifier_->DoOneRound(log_round_info)); RETURN_IF_FALSE(RemoveFixedAndEquivalentVariables(log_round_info)); + RETURN_IF_FALSE(LevelZeroPropagate()); // TODO(user): Add a small wrapper function to time this. RETURN_IF_FALSE(LevelZeroPropagate()); - sat_solver_->MinimizeSomeClauses(/*decisions_budget=*/1000); - RETURN_IF_FALSE(LevelZeroPropagate()); - - RETURN_IF_FALSE(SubsumeAndStrenghtenRound(log_round_info)); + const auto old_counter = sat_solver_->counters(); + RETURN_IF_FALSE(sat_solver_->MinimizeByPropagation()); + const int64_t mini_num_clause = + sat_solver_->counters().minimization_num_clauses - + old_counter.minimization_num_clauses; + const int64_t mini_num_removed = + sat_solver_->counters().minimization_num_removed_literals - + old_counter.minimization_num_removed_literals; RETURN_IF_FALSE(RemoveFixedAndEquivalentVariables(log_round_info)); - blocked_clause_simplifier_->DoOneRound(log_round_info); - RETURN_IF_FALSE(bounded_variable_elimination_->DoOneRound(log_round_info)); + RETURN_IF_FALSE(SubsumeAndStrenghtenRound(log_round_info)); + RETURN_IF_FALSE(RemoveFixedAndEquivalentVariables(log_round_info)); + + // TODO(user): try to enable these? The problem is that we can only remove + // variables not used the non-pure SAT part of a model. + if (/*DISABLES_CODE*/ (false)) { + blocked_clause_simplifier_->DoOneRound(log_round_info); + RETURN_IF_FALSE(bounded_variable_elimination_->DoOneRound(log_round_info)); + } RETURN_IF_FALSE(LevelZeroPropagate()); total_dtime_ += time_limit_->GetElapsedDeterministicTime() - start_dtime; - LOG_IF(INFO, log_info) - << "Presolve." - << " num_fixed: " << trail_->Index() - << " num_redundant: " << implication_graph_->num_redundant_literals() / 2 - << "/" << sat_solver_->NumVariables() - << " num_implications: " << implication_graph_->num_implications() - << " num_watched_clauses: " << clause_manager_->num_watched_clauses() - << " dtime: " << time_limit_->GetElapsedDeterministicTime() - start_dtime - << " wtime: " << wall_timer.Get() - << " non-probing time: " << (wall_timer.Get() - probing_time); + if (log_info) { + SOLVER_LOG( + logger_, "Inprocessing.", " fixed:", trail_->Index(), + " equiv:", implication_graph_->num_redundant_literals() / 2, + " bools:", sat_solver_->NumVariables(), + " implications:", implication_graph_->num_implications(), + " watched:", clause_manager_->num_watched_clauses(), + " minimization:", mini_num_clause, "|", mini_num_removed, + " dtime:", time_limit_->GetElapsedDeterministicTime() - start_dtime, + " wtime:", wall_timer.Get(), + " np_wtime:", (wall_timer.Get() - probing_time)); + } + DCHECK_EQ(sat_solver_->CurrentDecisionLevel(), 0); decision_policy_->MaybeEnablePhaseSaving(/*save_phase=*/true); return true; } @@ -399,8 +425,9 @@ bool Inprocessing::SubsumeAndStrenghtenRound(bool log_info) { // TODO(user): probably faster without the size indirection. std::vector clauses = clause_manager_->AllClausesInCreationOrder(); - std::sort(clauses.begin(), clauses.end(), - [](SatClause* a, SatClause* b) { return a->size() < b->size(); }); + std::stable_sort( + clauses.begin(), clauses.end(), + [](SatClause* a, SatClause* b) { return a->size() < b->size(); }); // Used to mark clause literals. const LiteralIndex num_literals(sat_solver_->NumVariables() * 2); @@ -516,16 +543,10 @@ bool Inprocessing::SubsumeAndStrenghtenRound(bool log_info) { if (!candidates_for_removal.empty()) { new_clause.clear(); for (const Literal l : clause->AsSpan()) { + if (l == candidates_for_removal[0]) continue; new_clause.push_back(l); } - - int new_size = 0; - for (const Literal l : new_clause) { - if (l == candidates_for_removal[0]) continue; - new_clause[new_size++] = l; - } - CHECK_EQ(new_size + 1, new_clause.size()); - new_clause.resize(new_size); + CHECK_EQ(new_clause.size() + 1, clause->size()); num_removed_literals += clause->size() - new_clause.size(); if (!clause_manager_->InprocessingRewriteClause(clause, new_clause)) { @@ -616,7 +637,6 @@ bool StampingSimplifier::DoOneRound(bool log_info) { // Note that num_removed_literals_ do not count the literals of the subsumed // clauses. time_limit_->AdvanceDeterministicTime(dtime_); - log_info |= VLOG_IS_ON(1); LOG_IF(INFO, log_info) << "Stamping. num_removed_literals: " << num_removed_literals_ << " num_subsumed: " << num_subsumed_clauses_ @@ -642,7 +662,6 @@ bool StampingSimplifier::ComputeStampsForNextRound(bool log_info) { // TODO(user): compute some dtime, it is always zero currently. time_limit_->AdvanceDeterministicTime(dtime_); - log_info |= VLOG_IS_ON(1); LOG_IF(INFO, log_info) << "Prestamping." << " num_fixed: " << num_fixed_ << " dtime: " << dtime_ << " wtime: " << wall_timer.Get(); @@ -804,7 +823,7 @@ bool StampingSimplifier::ProcessClauses() { entries.push_back({i, true, first_stamps_[span[i].NegatedIndex()], last_stamps_[span[i].NegatedIndex()]}); } - if (clause->empty()) continue; + if (clause->IsRemoved()) continue; // The sort should be dominant. if (!entries.empty()) { @@ -875,7 +894,7 @@ bool StampingSimplifier::ProcessClauses() { } } - if (clause->empty()) continue; + if (clause->IsRemoved()) continue; // Strengthen the clause. if (!to_remove.empty() || entries.size() < span.size()) { @@ -918,7 +937,10 @@ void BlockedClauseSimplifier::DoOneRound(bool log_info) { const Literal l = queue_.front(); in_queue_[l] = false; queue_.pop_front(); - ProcessLiteral(l); + + // Avoid doing too much work here on large problem. + // Note that we still what to empty the queue. + if (num_inspected_literals_ <= 1e9) ProcessLiteral(l); } // Release some memory. @@ -997,7 +1019,7 @@ void BlockedClauseSimplifier::ProcessLiteral(Literal current_literal) { // clauses_to_process clauses that resolve trivially with that clause. std::vector clauses_to_process; for (const ClauseIndex i : literal_to_clauses_[current_literal]) { - if (clauses_[i]->empty()) continue; + if (clauses_[i]->IsRemoved()) continue; // Blocked with respect to binary clause only? all marked binary should have // their negation in clause. @@ -1067,7 +1089,7 @@ bool BlockedClauseSimplifier::ClauseIsBlocked( // all clauses that are used in a non-blocked certificate first in the list. for (const ClauseIndex i : literal_to_clauses_[current_literal.NegatedIndex()]) { - if (clauses_[i]->empty()) continue; + if (clauses_[i]->IsRemoved()) continue; bool some_marked = false; for (const Literal l : clauses_[i]->AsSpan()) { // TODO(user): we can be faster here by only updating it at the end? @@ -1134,7 +1156,6 @@ bool BoundedVariableElimination::DoOneRound(bool log_info) { queue_.Reserve(num_variables); for (BooleanVariable v(0); v < num_variables; ++v) { if (assignment_.VariableIsAssigned(v)) continue; - if (implication_graph_->IsRemoved(Literal(v, true))) continue; UpdatePriorityQueue(v); } @@ -1173,13 +1194,12 @@ bool BoundedVariableElimination::DoOneRound(bool log_info) { need_to_be_updated_.clear(); } - implication_graph_->RemoveFixedVariables(); - implication_graph_->CleanupAllRemovedVariables(); + if (!Propagate()) return false; + implication_graph_->CleanupAllRemovedAndFixedVariables(); // Remove all redundant clause containing a removed literal. This avoid to // re-introduce a removed literal via conflict learning. for (SatClause* c : clause_manager_->AllClausesInCreationOrder()) { - if (!clause_manager_->IsRemovable(c)) continue; bool remove = false; for (const Literal l : c->AsSpan()) { if (implication_graph_->IsRemoved(l)) { @@ -1229,7 +1249,7 @@ bool BoundedVariableElimination::RemoveLiteralFromClause( if (!clause_manager_->InprocessingRewriteClause(sat_clause, resolvant_)) { return false; } - if (sat_clause->empty()) { + if (sat_clause->IsRemoved()) { --num_clauses_diff_; for (const Literal l : resolvant_) literal_to_num_clauses_[l]--; } else { @@ -1245,14 +1265,14 @@ bool BoundedVariableElimination::Propagate() { const Literal l = (*trail_)[propagation_index_]; for (const ClauseIndex index : literal_to_clauses_[l]) { - if (clauses_[index]->empty()) continue; + if (clauses_[index]->IsRemoved()) continue; num_clauses_diff_--; num_literals_diff_ -= clauses_[index]->size(); clause_manager_->InprocessingRemoveClause(clauses_[index]); } literal_to_clauses_[l].clear(); for (const ClauseIndex index : literal_to_clauses_[l.NegatedIndex()]) { - if (clauses_[index]->empty()) continue; + if (clauses_[index]->IsRemoved()) continue; if (!RemoveLiteralFromClause(l.Negated(), clauses_[index])) return false; } literal_to_clauses_[l.NegatedIndex()].clear(); @@ -1270,6 +1290,8 @@ int BoundedVariableElimination::NumClausesContaining(Literal l) { // TODO(user): Only enqueue variable that can be removed. void BoundedVariableElimination::UpdatePriorityQueue(BooleanVariable var) { if (assignment_.VariableIsAssigned(var)) return; + if (implication_graph_->IsRemoved(Literal(var, true))) return; + if (implication_graph_->IsRedundant(Literal(var, true))) return; const int priority = -NumClausesContaining(Literal(var, true)) - NumClausesContaining(Literal(var, false)); if (queue_.Contains(var.value())) { @@ -1496,6 +1518,7 @@ bool BoundedVariableElimination::CrossProduct(BooleanVariable var) { const Literal lit(var, true); const Literal not_lit(var, false); + DCHECK(!implication_graph_->IsRedundant(lit)); { const int s1 = NumClausesContaining(lit); const int s2 = NumClausesContaining(not_lit); @@ -1512,15 +1535,6 @@ bool BoundedVariableElimination::CrossProduct(BooleanVariable var) { DeleteAllClausesContaining(not_lit); return true; } - if (implication_graph_->IsRedundant(lit)) { - // TODO(user): do that elsewhere? - CHECK_EQ(s1, 1); - CHECK_EQ(s2, 1); - CHECK_EQ(implication_graph_->NumImplicationOnVariableRemoval(var), 0); - num_eliminated_variables_++; - implication_graph_->RemoveBooleanVariable(var, &postsolve_->clauses); - return true; - } // Heuristic. Abort if the work required to decide if var should be removed // seems to big. diff --git a/ortools/sat/sat_inprocessing.h b/ortools/sat/sat_inprocessing.h index c5a44c9bd0..36f0c2e703 100644 --- a/ortools/sat/sat_inprocessing.h +++ b/ortools/sat/sat_inprocessing.h @@ -108,6 +108,7 @@ class Inprocessing { model->GetOrCreate()), bounded_variable_elimination_( model->GetOrCreate()), + postsolve_(model->GetOrCreate()), logger_(model->GetOrCreate()), model_(model) {} @@ -155,6 +156,7 @@ class Inprocessing { StampingSimplifier* stamping_simplifier_; BlockedClauseSimplifier* blocked_clause_simplifier_; BoundedVariableElimination* bounded_variable_elimination_; + PostsolveClauses* postsolve_; SolverLogger* logger_; double total_dtime_ = 0.0; diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 68e0bd1a15..5c74e20e4a 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -23,7 +23,7 @@ option csharp_namespace = "Google.OrTools.Sat"; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 269 +// NEXT TAG: 272 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -194,18 +194,21 @@ message SatParameters { // learnt clause minimization approach for CDCL Sat Solvers", // https://www.ijcai.org/proceedings/2017/0098.pdf // - // For now, we have a somewhat simpler implementation where every x restart we - // spend y decisions on clause minimization. The minimization technique is the - // same as the one used to minimize core in max-sat. We also minimize problem - // clauses and not just the learned clause that we keep forever like in the - // paper. + // For now, we have a somewhat simpler implementation where, each time we + // are at level 0, if more that x num decision where taken, we + // spend some y decision on clause minimization. We try to keep y / x equal + // to the given ratio. + // + // The minimization technique is the same as the one used to minimize core in + // max-sat. We also minimize problem clauses and not just the learned clause + // that we keep forever like in the paper. // // Changing these parameters or the kind of clause we minimize seems to have // a big impact on the overall perf on our benchmarks. So this technique seems // definitely useful, but it is hard to tune properly. // - // A negative period disable the feature. - optional int32 minimize_with_propagation_restart_period = 96 [default = -1]; + // A ratio of zero disable the feature. + optional double minimize_with_propagation_ratio = 96 [default = 0.0]; optional int32 minimize_with_propagation_num_decisions = 97 [default = 1000]; // ========================================================================== @@ -811,6 +814,12 @@ message SatParameters { optional bool use_energetic_reasoning_in_no_overlap_2d = 213 [default = false]; + // When this is true, the no_overlap_2d constraint is reinforced with + // an energetic reasoning that uses an area-based energy. This can be combined + // with the two other overlap heuristics above. + optional bool use_area_energetic_reasoning_in_no_overlap_2d = 271 + [default = false]; + // Performs an extra step of propagation in the no_overlap_2d constraint by // looking at all pairs of intervals. optional bool use_pairwise_reasoning_in_no_overlap_2d = 251 [default = false]; @@ -1027,16 +1036,14 @@ message SatParameters { // branch on the value that lead to the best objective first. optional bool exploit_objective = 131 [default = true]; - // If set at zero (the default), it is disabled. Otherwise the solver attempts - // probing at every 'probing_period' root node. Period of 1 enables probing at - // every root node. - optional int64 probing_period_at_root = 142 [default = 0]; - // If true, search will continuously probe Boolean variables, and integer // variable bounds. This parameter is set to true in parallel on the probing // worker. optional bool use_probing_search = 176 [default = false]; + // Use extended probing (probe bool_or, at_most_one, exactly_one). + optional bool use_extended_probing = 269 [default = true]; + // Add a shaving phase (where the solver tries to prove that the lower or // upper bound of a variable are infeasible) to the probing search. optional bool use_shaving_in_probing_search = 204 [default = true]; @@ -1321,6 +1328,10 @@ message SatParameters { // intervals, but 1M intervals in the no-overlap constraints covering them. optional bool use_combined_no_overlap = 133 [default = false]; + // All at_most_one constraints with a size <= param will be replaced by a + // quadratic number of binary implications. + optional int32 at_most_one_max_expansion_size = 270 [default = 9]; + // Indicates if the CP-SAT layer should catch Control-C (SIGINT) signals // when calling solve. If set, catching the SIGINT signal will terminate the // search gracefully, as if a time limit was reached. diff --git a/ortools/sat/sat_solver.cc b/ortools/sat/sat_solver.cc index c7f7a85149..be781e1b35 100644 --- a/ortools/sat/sat_solver.cc +++ b/ortools/sat/sat_solver.cc @@ -1006,7 +1006,7 @@ SatSolver::Status SatSolver::EnqueueDecisionAndBacktrackOnConflict( bool SatSolver::EnqueueDecisionIfNotConflicting(Literal true_literal) { SCOPED_TIME_STAT(&stats_); - CHECK(PropagationIsDone()); + DCHECK(PropagationIsDone()); if (model_is_unsat_) return kUnsatTrailIndex; const int current_level = CurrentDecisionLevel(); @@ -1128,6 +1128,9 @@ void SatSolver::KeepAllClauseUsedToInfer(BooleanVariable variable) { // Keep this clause. clauses_propagator_->mutable_clauses_info()->erase(clause); } + if (trail_->AssignmentType(var) == AssignmentType::kSearchDecision) { + continue; + } for (const Literal l : trail_->Reason(var)) { const AssignmentInfo& info = trail_->Info(l.Variable()); if (info.level == 0) continue; @@ -1160,10 +1163,8 @@ bool SatSolver::SubsumptionIsInteresting(BooleanVariable variable) { const BooleanVariable var = (*trail_)[trail_index].Variable(); const int type = trail_->AssignmentType(var); - if (type != AssignmentType::kSearchDecision && type != binary_id && - type != clause_id) { - return false; - } + if (type == AssignmentType::kSearchDecision) continue; + if (type != binary_id && type != clause_id) return false; SatClause* clause = ReasonClauseOrNull(var); if (clause != nullptr && clauses_propagator_->IsRemovable(clause)) { ++num_clause_to_mark_as_non_deletable; @@ -1185,10 +1186,12 @@ bool SatSolver::SubsumptionIsInteresting(BooleanVariable variable) { // Ideally this should be scheduled after other faster in-processing technique. void SatSolver::TryToMinimizeClause(SatClause* clause) { CHECK_EQ(CurrentDecisionLevel(), 0); + CHECK(clause != nullptr); ++counters_.minimization_num_clauses; absl::btree_set moved_last; std::vector candidate(clause->begin(), clause->end()); + while (!model_is_unsat_) { // We want each literal in candidate to appear last once in our propagation // order. We want to do that while maximizing the reutilization of the @@ -1241,7 +1244,7 @@ void SatSolver::TryToMinimizeClause(SatClause* clause) { } else { ++counters_.minimization_num_decisions; EnqueueDecisionAndBackjumpOnConflict(literal.Negated()); - if (!clause->IsAttached()) { + if (clause->IsRemoved()) { Backtrack(0); return; } @@ -1321,12 +1324,7 @@ SatSolver::Status SatSolver::SolveInternal(TimeLimit* time_limit, } // Used to trigger clause minimization via propagation. - int64_t next_minimization_num_restart = - restart_->NumRestarts() + - parameters_->minimize_with_propagation_restart_period(); - if (parameters_->minimize_with_propagation_restart_period() < 0) { - next_minimization_num_restart = std::numeric_limits::max(); - } + ResetMinimizationByPropagationThreshold(); // Variables used to show the search progress. const int64_t kDisplayFrequency = 10000; @@ -1404,20 +1402,8 @@ SatSolver::Status SatSolver::SolveInternal(TimeLimit* time_limit, } // Clause minimization using propagation. - if (CurrentDecisionLevel() == 0 && - restart_->NumRestarts() >= next_minimization_num_restart) { - next_minimization_num_restart = - restart_->NumRestarts() + - parameters_->minimize_with_propagation_restart_period(); - MinimizeSomeClauses( - parameters_->minimize_with_propagation_num_decisions()); - - // Corner case: the minimization above being based on propagation may - // fix the remaining variables or prove UNSAT. - if (model_is_unsat_) return StatusWithLog(INFEASIBLE); - if (trail_->Index() == num_variables_.value()) { - return StatusWithLog(FEASIBLE); - } + if (CurrentDecisionLevel() == 0 && !MaybeMinimizeByPropagation()) { + return StatusWithLog(UnsatStatus()); } DCHECK_GE(CurrentDecisionLevel(), assumption_level_); @@ -1426,18 +1412,42 @@ SatSolver::Status SatSolver::SolveInternal(TimeLimit* time_limit, } } -void SatSolver::MinimizeSomeClauses(int decisions_budget) { +void SatSolver::ResetMinimizationByPropagationThreshold() { + if (parameters_->minimize_with_propagation_ratio() <= 0.0 || + parameters_->minimize_with_propagation_num_decisions() <= 0) { + minimization_by_propagation_threshold_ = + std::numeric_limits::max(); + return; + } + minimization_by_propagation_threshold_ = + counters_.num_branches + + static_cast( + static_cast( + parameters_->minimize_with_propagation_num_decisions()) / + parameters_->minimize_with_propagation_ratio()); +} + +bool SatSolver::MaybeMinimizeByPropagation() { + if (counters_.num_branches <= minimization_by_propagation_threshold_) { + return true; + } + return MinimizeByPropagation(); +} + +bool SatSolver::MinimizeByPropagation() { // Tricky: we don't want TryToMinimizeClause() to delete to_minimize // while we are processing it. block_clause_deletion_ = true; - const int64_t target_num_branches = counters_.num_branches + decisions_budget; + const int64_t target_num_branches = + counters_.num_branches + + parameters_->minimize_with_propagation_num_decisions(); while (counters_.num_branches < target_num_branches && (time_limit_ == nullptr || !time_limit_->LimitReached())) { SatClause* to_minimize = clauses_propagator_->NextClauseToMinimize(); if (to_minimize != nullptr) { TryToMinimizeClause(to_minimize); - if (model_is_unsat_) return; + if (model_is_unsat_) return false; } else { if (to_minimize == nullptr) { VLOG(1) << "Minimized all clauses, restarting from first one."; @@ -1449,6 +1459,13 @@ void SatSolver::MinimizeSomeClauses(int decisions_budget) { block_clause_deletion_ = false; clauses_propagator_->DeleteRemovedClauses(); + + ResetMinimizationByPropagationThreshold(); + + // Note(user): In some corner cases, the function above might find a + // feasible assignment. I think it is okay to ignore this special case + // that should only happen on trivial problems and just reset the solver. + return ResetToLevelZero(); } std::vector SatSolver::GetLastIncompatibleDecisions() { @@ -1747,7 +1764,7 @@ void SatSolver::ProcessNewlyFixedVariables() { // others. Note that none of the clause should be all false because we should // have detected a conflict before this is called. for (SatClause* clause : clauses_propagator_->AllClausesInCreationOrder()) { - if (!clause->IsAttached()) continue; + if (clause->IsRemoved()) continue; const size_t old_size = clause->size(); if (clause->RemoveFixedLiteralsAndTestIfTrue(trail_->Assignment())) { diff --git a/ortools/sat/sat_solver.h b/ortools/sat/sat_solver.h index eb4b211287..2b83ffa942 100644 --- a/ortools/sat/sat_solver.h +++ b/ortools/sat/sat_solver.h @@ -28,6 +28,7 @@ #include #include +#include "absl/base/attributes.h" #include "absl/container/flat_hash_map.h" #include "absl/log/check.h" #include "absl/strings/string_view.h" @@ -473,9 +474,13 @@ class SatSolver { // Mainly visible for testing. ABSL_MUST_USE_RESULT bool Propagate(); - // This must be called at level zero. It will spend the given num decision and - // use propagation to try to minimize some clauses from the database. - void MinimizeSomeClauses(int decisions_budget); + // This must be called at level zero. If enough decision were taken since the + // last ResetMinimizationByPropagationThreshold(), it will spend some decision + // (according to the parameters) and use propagation to try to minimize some + // clauses from the database. Return false on UNSAT. + bool MaybeMinimizeByPropagation(); + void ResetMinimizationByPropagationThreshold(); + bool MinimizeByPropagation(); // Sets the export function to the shared clauses manager. void SetShareBinaryClauseCallback(const std::function& @@ -684,7 +689,7 @@ class SatSolver { // function choose which one of the other functions to call depending on the // parameters. // - // Precondidtion: is_marked_ should be set to true for all the variables of + // Precondition: is_marked_ should be set to true for all the variables of // the conflict. It can also contains false non-conflict variables that // are implied by the negation of the 1-UIP conflict literal. void MinimizeConflict( @@ -701,7 +706,7 @@ class SatSolver { // - There is an unique literal with the highest decision level. // - This literal appears in the first position. // - All the other literals are of smaller decision level. - // - Ther is no literal with a decision level of zero. + // - There is no literal with a decision level of zero. bool IsConflictValid(const std::vector& literals); // Given the learned clause after a conflict, this computes the correct @@ -828,6 +833,8 @@ class SatSolver { // deleted. When it reaches zero, a clause cleanup is triggered. int num_learned_clause_before_cleanup_ = 0; + int64_t minimization_by_propagation_threshold_ = 0; + // Temporary members used during conflict analysis. SparseBitset is_marked_; SparseBitset is_independent_; diff --git a/ortools/sat/simplification.cc b/ortools/sat/simplification.cc index 1acc4d4aff..e1516db3c8 100644 --- a/ortools/sat/simplification.cc +++ b/ortools/sat/simplification.cc @@ -305,10 +305,10 @@ bool SatPresolver::ProcessAllClauses() { // Because on large problem we don't have a budget to process all clauses, // lets start by the smallest ones first. - std::sort(clause_to_process_.begin(), clause_to_process_.end(), - [this](ClauseIndex c1, ClauseIndex c2) { - return clauses_[c1].size() < clauses_[c2].size(); - }); + std::stable_sort(clause_to_process_.begin(), clause_to_process_.end(), + [this](ClauseIndex c1, ClauseIndex c2) { + return clauses_[c1].size() < clauses_[c2].size(); + }); while (!clause_to_process_.empty()) { const ClauseIndex ci = clause_to_process_.front(); in_clause_to_process_[ci] = false; @@ -381,6 +381,7 @@ bool SatPresolver::Presolve(const std::vector& can_be_removed) { return true; } +// TODO(user): Put work limit in place ! void SatPresolver::PresolveWithBva() { var_pq_elements_.clear(); // so we don't update it. InitializeBvaPriorityQueue(); @@ -388,6 +389,7 @@ void SatPresolver::PresolveWithBva() { const LiteralIndex lit = bva_pq_.Top()->literal; bva_pq_.Pop(); SimpleBva(lit); + if (time_limit_ != nullptr && time_limit_->LimitReached()) break; } } @@ -403,7 +405,7 @@ void SatPresolver::SimpleBva(LiteralIndex l) { m_cls_ = literal_to_clauses_[l]; int reduction = 0; - while (true) { + for (int loop = 0; loop < 100; ++loop) { LiteralIndex lmax = kNoLiteralIndex; int max_size = 0;