diff --git a/ortools/sat/2d_mandatory_overlap_propagator.cc b/ortools/sat/2d_mandatory_overlap_propagator.cc index 9738f86507..6974d67f35 100644 --- a/ortools/sat/2d_mandatory_overlap_propagator.cc +++ b/ortools/sat/2d_mandatory_overlap_propagator.cc @@ -24,6 +24,7 @@ #include "absl/types/span.h" #include "ortools/sat/diffn_util.h" #include "ortools/sat/integer.h" +#include "ortools/sat/integer_base.h" #include "ortools/sat/model.h" #include "ortools/sat/no_overlap_2d_helper.h" #include "ortools/sat/scheduling_helpers.h" @@ -57,18 +58,22 @@ bool MandatoryOverlapPropagator::Propagate() { bool has_zero_area_boxes = false; absl::Span tasks = helper_.x_helper().TaskByIncreasingNegatedStartMax(); - for (int i = tasks.size() - 1; i >= 0; --i) { + for (int i = tasks.size(); --i >= 0;) { const int b = tasks[i].task_index; if (!helper_.IsPresent(b)) continue; - const ItemWithVariableSize item = helper_.GetItemWithVariableSize(b); - if (item.x.start_max > item.x.end_min || - item.y.start_max > item.y.end_min) { - continue; - } - mandatory_regions_.push_back({.x_min = item.x.start_max, - .x_max = item.x.end_min, - .y_min = item.y.start_max, - .y_max = item.y.end_min}); + + const IntegerValue x_start_max = helper_.x_helper().StartMax(b); + const IntegerValue x_end_min = helper_.x_helper().EndMin(b); + if (x_start_max > x_end_min) continue; + + const IntegerValue y_start_max = helper_.y_helper().StartMax(b); + const IntegerValue y_end_min = helper_.y_helper().EndMin(b); + if (y_start_max > y_end_min) continue; + + mandatory_regions_.push_back({.x_min = x_start_max, + .x_max = x_end_min, + .y_min = y_start_max, + .y_max = y_end_min}); mandatory_regions_index_.push_back(b); if (mandatory_regions_.back().SizeX() == 0 || diff --git a/ortools/sat/2d_try_edge_propagator.cc b/ortools/sat/2d_try_edge_propagator.cc index 44ba155d2d..c8616d9fda 100644 --- a/ortools/sat/2d_try_edge_propagator.cc +++ b/ortools/sat/2d_try_edge_propagator.cc @@ -23,6 +23,7 @@ #include "absl/log/check.h" #include "absl/log/log.h" #include "absl/log/vlog_is_on.h" +#include "absl/types/span.h" #include "ortools/base/stl_util.h" #include "ortools/sat/diffn_util.h" #include "ortools/sat/integer.h" @@ -31,6 +32,7 @@ #include "ortools/sat/no_overlap_2d_helper.h" #include "ortools/sat/synchronization.h" #include "ortools/sat/util.h" +#include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_heuristics.h" #include "ortools/set_cover/set_cover_invariant.h" #include "ortools/set_cover/set_cover_model.h" @@ -117,9 +119,10 @@ bool TryEdgeRectanglePropagator::CanPlace( .x_max = position.first + active_box_ranges_[box_index].x_size, .y_min = position.second, .y_max = position.second + active_box_ranges_[box_index].y_size}; + const absl::Span mandatory_regions = mandatory_regions_; for (const int i : has_mandatory_region_) { if (i == box_index) continue; - const Rectangle& mandatory_region = mandatory_regions_[i]; + const Rectangle& mandatory_region = mandatory_regions[i]; if (!mandatory_region.IsDisjoint(placed_box)) { if (with_conflict != nullptr) { with_conflict->AppendToLastVector(i); @@ -149,10 +152,12 @@ bool TryEdgeRectanglePropagator::Propagate() { // If a mandatory region is changed, we need to replace any cached box that // now became overlapping with it. + const int num_active_ranges = active_box_ranges_.size(); for (const int mandatory_idx : changed_mandatory_) { - for (int i = 0; i < active_box_ranges_.size(); i++) { + const Rectangle& mandatory_region = mandatory_regions_[mandatory_idx]; + for (int i = 0; i < num_active_ranges; i++) { if (i == mandatory_idx || !is_in_cache_[i]) continue; - if (!placed_boxes_[i].IsDisjoint(mandatory_regions_[mandatory_idx])) { + if (!placed_boxes_[i].IsDisjoint(mandatory_region)) { changed_item_.push_back(i); is_in_cache_[i] = false; } @@ -323,8 +328,10 @@ std::vector TryEdgeRectanglePropagator::GetMinimumProblemWithPropagation( } DCHECK(model.ComputeFeasibility()); SetCoverInvariant inv(&model); - GreedySolutionGenerator greedy_search(&inv); + LazyElementDegreeSolutionGenerator greedy_search(&inv); CHECK(greedy_search.NextSolution()); + LazySteepestSearch steepest_search(&inv); + CHECK(steepest_search.NextSolution()); GuidedLocalSearch search(&inv); CHECK(search.SetMaxIterations(100).NextSolution()); DCHECK(inv.CheckConsistency( diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index 4d80a69db6..eba29fa24c 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -101,6 +101,7 @@ cc_library( deps = [ ":diffn_util", ":integer", + ":integer_base", ":model", ":no_overlap_2d_helper", ":scheduling_helpers", @@ -2377,19 +2378,23 @@ cc_test( size = "small", srcs = [ "boolean_problem_test.cc", - "opb_reader.h", ], deps = [ ":boolean_problem", - ":boolean_problem_cc_proto", + ":cp_model_cc_proto", + ":cp_model_checker", + ":cp_model_symmetries", + ":opb_reader", + ":sat_parameters_cc_proto", "//ortools/algorithms:sparse_permutation", "//ortools/base", "//ortools/base:file", "//ortools/base:gmock_main", "//ortools/base:path", "//ortools/util:filelineiter", + "//ortools/util:logging", + "//ortools/util:time_limit", "@abseil-cpp//absl/log:check", - "@abseil-cpp//absl/status", "@abseil-cpp//absl/strings", ], ) @@ -3320,6 +3325,7 @@ cc_library( ":synchronization", ":util", "//ortools/base:stl_util", + "//ortools/set_cover:base_types", "//ortools/set_cover:set_cover_heuristics", "//ortools/set_cover:set_cover_invariant", "//ortools/set_cover:set_cover_model", @@ -3328,6 +3334,7 @@ cc_library( "@abseil-cpp//absl/log", "@abseil-cpp//absl/log:check", "@abseil-cpp//absl/log:vlog_is_on", + "@abseil-cpp//absl/types:span", ], ) @@ -3737,19 +3744,35 @@ cc_library( ], ) +cc_library( + name = "opb_reader", + hdrs = ["opb_reader.h"], + deps = [ + ":cp_model_cc_proto", + ":cp_model_utils", + "//ortools/base", + "//ortools/base:stl_util", + "//ortools/util:filelineiter", + "@abseil-cpp//absl/container:flat_hash_map", + "@abseil-cpp//absl/container:flat_hash_set", + "@abseil-cpp//absl/log", + "@abseil-cpp//absl/log:check", + "@abseil-cpp//absl/strings", + "@abseil-cpp//absl/types:span", + ], +) + cc_binary( name = "sat_runner", srcs = [ - "opb_reader.h", "sat_runner.cc", ], deps = [ - ":boolean_problem", - ":boolean_problem_cc_proto", ":cp_model_cc_proto", ":cp_model_solver", ":cp_model_utils", ":model", + ":opb_reader", ":sat_cnf_reader", ":sat_parameters_cc_proto", "//ortools/base", diff --git a/ortools/sat/boolean_problem_test.cc b/ortools/sat/boolean_problem_test.cc index 2f086c2381..9d68a9e999 100644 --- a/ortools/sat/boolean_problem_test.cc +++ b/ortools/sat/boolean_problem_test.cc @@ -18,7 +18,6 @@ #include #include "absl/log/check.h" -#include "absl/status/status.h" #include "absl/strings/str_cat.h" #include "absl/strings/string_view.h" #include "gtest/gtest.h" @@ -26,66 +25,97 @@ #include "ortools/base/helpers.h" #include "ortools/base/options.h" #include "ortools/base/path.h" -#include "ortools/sat/boolean_problem.pb.h" +#include "ortools/sat/cp_model.pb.h" +#include "ortools/sat/cp_model_checker.h" +#include "ortools/sat/cp_model_symmetries.h" #include "ortools/sat/opb_reader.h" +#include "ortools/sat/sat_parameters.pb.h" +#include "ortools/util/logging.h" +#include "ortools/util/time_limit.h" namespace operations_research { namespace sat { namespace { -TEST(ValidateBooleanProblemTest, Ok) { +TEST(LoadAndValidateBooleanProblemTest, Ok) { std::string file = "min: 1 x1 1 x2 ;\n" "1 x1 1 x2 >= 1 ;\n" "1 x1 1 x2 >= 1 ;\n"; - LinearBooleanProblem problem; + CpModelProto problem; OpbReader reader; const std::string filename = file::JoinPath(::testing::TempDir(), "file.opb"); CHECK_OK(file::SetContents(filename, file, file::Defaults())); - CHECK(reader.Load(filename, &problem)); - EXPECT_TRUE(ValidateBooleanProblem(problem).ok()); + CHECK(reader.LoadAndValidate(filename, &problem)); + EXPECT_TRUE(ValidateCpModel(problem).empty()); + EXPECT_EQ(problem.variables_size(), 2); + EXPECT_EQ(reader.num_variables(), 2); } -TEST(ValidateBooleanProblemTest, ZeroCoefficients) { +TEST(LoadAndValidateBooleanProblemTest, NonLinear) { + std::string file = + "min: 1 x1 1 x2 ~x3;\n" + "1 x1 1 x2 >= 1 ;\n" + "1 x1 x2 x3 >= 1 ;\n"; + CpModelProto problem; + OpbReader reader; + const std::string filename = file::JoinPath(::testing::TempDir(), "file.opb"); + CHECK_OK(file::SetContents(filename, file, file::Defaults())); + CHECK(reader.LoadAndValidate(filename, &problem)); + EXPECT_TRUE(ValidateCpModel(problem).empty()); + EXPECT_EQ(problem.variables_size(), 5); + EXPECT_EQ(reader.num_variables(), 3); +} + +TEST(LoadAndValidateBooleanProblemTest, Weighted) { + std::string file = + "1 x1 1 x2 >= 1 ;\n" + "[32] 2 x1 = 2 ;\n"; + CpModelProto problem; + OpbReader reader; + const std::string filename = file::JoinPath(::testing::TempDir(), "file.opb"); + CHECK_OK(file::SetContents(filename, file, file::Defaults())); + CHECK(reader.LoadAndValidate(filename, &problem)); + EXPECT_TRUE(ValidateCpModel(problem).empty()); + EXPECT_EQ(problem.variables_size(), 3); + EXPECT_EQ(reader.num_variables(), 2); + ASSERT_TRUE(problem.has_objective()); + EXPECT_EQ(problem.objective().vars_size(), 1); + EXPECT_EQ(problem.objective().vars(0), 2); + EXPECT_EQ(problem.objective().coeffs(0), 32); +} + +TEST(LoadAndValidateBooleanProblemTest, ZeroCoefficients) { std::string file = "min: 1 x1 1 x2 ;\n" "1 x1 0 x2 >= 1 ;\n" "1 x1 1 x2 >= 1 ;\n"; - LinearBooleanProblem problem; + CpModelProto problem; OpbReader reader; const std::string filename = file::JoinPath(::testing::TempDir(), "file2.opb"); CHECK_OK(file::SetContents(filename, file, file::Defaults())); - CHECK(reader.Load(filename, &problem)); - EXPECT_FALSE(ValidateBooleanProblem(problem).ok()); -} - -TEST(ValidateBooleanProblemTest, DuplicateEntries) { - std::string file = - "min: 1 x1 1 x2 ;\n" - "1 x1 1 x2 1 x1 >= 1 ;\n" - "1 x1 1 x2 >= 1 ;\n"; - LinearBooleanProblem problem; - OpbReader reader; - const std::string filename = - file::JoinPath(::testing::TempDir(), "file3.opb"); - CHECK_OK(file::SetContents(filename, file, file::Defaults())); - CHECK(reader.Load(filename, &problem)); - EXPECT_FALSE(ValidateBooleanProblem(problem).ok()); + EXPECT_FALSE(reader.LoadAndValidate(filename, &problem)); } void FindSymmetries( absl::string_view file, std::vector>* generators) { - LinearBooleanProblem problem; + CpModelProto problem; OpbReader reader; static int counter = 4; ++counter; const std::string filename = file::JoinPath( ::testing::TempDir(), absl::StrCat("file", counter, ".opb")); CHECK_OK(file::SetContents(filename, file, file::Defaults())); - CHECK(reader.Load(filename, &problem)); - FindLinearBooleanProblemSymmetries(problem, generators); + CHECK(reader.LoadAndValidate(filename, &problem)); + SatParameters params; + params.set_symmetry_level(2); + params.set_symmetry_detection_deterministic_time_limit(1000); + SolverLogger logger; + logger.EnableLogging(true); + TimeLimit time_limit; + FindCpModelSymmetries(params, problem, generators, &logger, &time_limit); } TEST(FindLinearBooleanProblemSymmetriesTest, ProblemWithSymmetry1) { @@ -96,19 +126,17 @@ TEST(FindLinearBooleanProblemSymmetriesTest, ProblemWithSymmetry1) { std::vector> generators; FindSymmetries(file, &generators); - // Note that the permutation is on the literals: - // xi maps to 2i and not(xi) maps to 2i + 1. EXPECT_EQ(generators.size(), 1); - EXPECT_EQ(generators[0]->DebugString(), "(0 2) (1 3)"); + EXPECT_EQ(generators[0]->DebugString(), "(0 1)"); } -TEST(FindLinearBooleanProblemSymmetriesTest, ProblemWithSymmetry2) { +TEST(DISABLED_FindLinearBooleanProblemSymmetriesTest, ProblemWithSymmetry2) { std::string file = "min: 1 x1 1 x2 ;\n" "-3 x1 -2 x2 >= -1 ;\n"; // This is simplified to both x1 and x2 false. std::vector> generators; FindSymmetries(file, &generators); - EXPECT_EQ(generators.size(), 1); + ASSERT_EQ(generators.size(), 1); EXPECT_EQ(generators[0]->DebugString(), "(0 2) (1 3)"); } @@ -120,18 +148,18 @@ TEST(FindLinearBooleanProblemSymmetriesTest, ProblemWithSymmetry3) { " 1 x3 2 x1 3 x2 >= 2 ;\n"; std::vector> generators; FindSymmetries(file, &generators); - EXPECT_EQ(generators.size(), 1); - EXPECT_EQ(generators[0]->DebugString(), "(0 4 2) (1 5 3)"); + ASSERT_EQ(generators.size(), 1); + EXPECT_EQ(generators[0]->DebugString(), "(0 2 1)"); } -TEST(FindLinearBooleanProblemSymmetriesTest, ProblemWithSymmetry4) { +TEST(DISABLED_FindLinearBooleanProblemSymmetriesTest, ProblemWithSymmetry4) { std::string file = "min: 1 x1;\n" " 1 x1 2 x2 >= 2 ;\n" " 1 x1 -2 x3 >= 0 ;\n"; std::vector> generators; FindSymmetries(file, &generators); - EXPECT_EQ(generators.size(), 1); + ASSERT_EQ(generators.size(), 1); // x2 and not(x3) are equivalent. EXPECT_EQ(generators[0]->DebugString(), "(2 5) (3 4)"); diff --git a/ortools/sat/clause.cc b/ortools/sat/clause.cc index 83c6ff26ae..217d949a58 100644 --- a/ortools/sat/clause.cc +++ b/ortools/sat/clause.cc @@ -574,7 +574,8 @@ bool BinaryImplicationGraph::AddBinaryClause(Literal a, Literal b) { // Tricky: If this is the first clause, the propagator will be added and // assumed to be in a "propagated" state. This makes sure this is the case. - if (IsEmpty()) propagation_trail_index_ = trail_->Index(); + if (no_constraint_ever_added_) propagation_trail_index_ = trail_->Index(); + no_constraint_ever_added_ = false; if (drat_proof_handler_ != nullptr) { // TODO(user): Like this we will duplicate all binary clause from the @@ -600,7 +601,6 @@ bool BinaryImplicationGraph::AddBinaryClause(Literal a, Literal b) { NotifyPossibleDuplicate(a); NotifyPossibleDuplicate(b); is_dag_ = false; - num_implications_ += 2; if (enable_sharing_ && add_binary_callback_ != nullptr) { add_binary_callback_(a, b); @@ -634,6 +634,10 @@ bool BinaryImplicationGraph::AddAtMostOne( DCHECK_EQ(trail_->CurrentDecisionLevel(), 0); if (at_most_one.size() <= 1) return true; + // Same as for AddBinaryClause(). + if (no_constraint_ever_added_) propagation_trail_index_ = trail_->Index(); + no_constraint_ever_added_ = false; + // Temporarily copy the at_most_one constraint at the end of // at_most_one_buffer_. It will be cleaned up and added by // CleanUpAndAddAtMostOnes(). @@ -788,7 +792,6 @@ bool BinaryImplicationGraph::CleanUpAndAddAtMostOnes(int base_index) { NotifyPossibleDuplicate(a); } } - num_implications_ += at_most_one.size() * (at_most_one.size() - 1); // This will erase the at_most_one from the buffer. local_end = local_start; @@ -1173,7 +1176,7 @@ void BinaryImplicationGraph::RemoveFixedVariables() { // TODO(user): This might be a bit slow. Do not call all the time if needed, // this shouldn't change the correctness of the code. at_most_ones_.clear(); - CleanUpAndAddAtMostOnes(/*base_index=*/0); + CHECK(CleanUpAndAddAtMostOnes(/*base_index=*/0)); DCHECK(InvariantsAreOk()); } @@ -1428,6 +1431,7 @@ bool BinaryImplicationGraph::DetectEquivalences(bool log_info) { const Literal rep = RepresentativeOf(l); if (rep.Index() != representative) representative_list.push_back(rep); } + dtime += 1e-8 * static_cast(ref.size()); // Add representative <=> literal. // @@ -1449,7 +1453,8 @@ bool BinaryImplicationGraph::DetectEquivalences(bool log_info) { // one. at_most_ones_.clear(); int saved_trail_index = propagation_trail_index_; - CleanUpAndAddAtMostOnes(/*base_index=*/0); + if (!CleanUpAndAddAtMostOnes(/*base_index=*/0)) return false; + // This might have run the propagation on a few variables without taking // into account the AMOs. Propagate again. // @@ -1459,12 +1464,6 @@ bool BinaryImplicationGraph::DetectEquivalences(bool log_info) { propagation_trail_index_ = saved_trail_index; Propagate(trail_); } - - num_implications_ = 0; - for (LiteralIndex i(0); i < size; ++i) { - num_implications_ += implications_[i].size(); - } - dtime += 2e-8 * num_implications_; } time_limit_->AdvanceDeterministicTime(dtime); @@ -1475,8 +1474,9 @@ bool BinaryImplicationGraph::DetectEquivalences(bool log_info) { LOG_IF(INFO, log_info) << "SCC. " << num_equivalences << " redundant equivalent literals. " << num_fixed_during_scc << " fixed. " - << num_implications_ << " implications left. " - << implications_.size() << " literals." + << ComputeNumImplicationsForLog() + << " implications left. " << implications_.size() + << " literals." << " size of at_most_one buffer = " << at_most_one_buffer_.size() << "." << " dtime: " << dtime @@ -1640,7 +1640,6 @@ bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { direct_implications.resize(new_size); direct_implications.shrink_to_fit(); num_new_redundant_implications += diff; - num_implications_ -= diff; // Abort if the computation involved is too big. if (work_done_in_mark_descendants_ > 1e8) { @@ -1681,7 +1680,8 @@ bool BinaryImplicationGraph::ComputeTransitiveReduction(bool log_info) { num_redundant_implications_ += num_new_redundant_implications; LOG_IF(INFO, log_info) << "Transitive reduction removed " << num_new_redundant_implications << " literals. " - << num_fixed << " fixed. " << num_implications_ + << num_fixed << " fixed. " + << ComputeNumImplicationsForLog() << " implications left. " << implications_.size() << " literals." << " dtime: " << dtime << " wtime: " << wall_timer.Get() @@ -2681,7 +2681,7 @@ void BinaryImplicationGraph::CleanupAllRemovedAndFixedVariables() { // Clean-up at most ones. at_most_ones_.clear(); - CleanUpAndAddAtMostOnes(/*base_index=*/0); + CHECK(CleanUpAndAddAtMostOnes(/*base_index=*/0)); // Note that to please the invariant() we also removed fixed literal above. DCHECK(InvariantsAreOk()); diff --git a/ortools/sat/clause.h b/ortools/sat/clause.h index 752acfff63..de837c5a26 100644 --- a/ortools/sat/clause.h +++ b/ortools/sat/clause.h @@ -520,10 +520,12 @@ class BinaryImplicationGraph : public SatPropagator { // Resizes the data structure. void Resize(int num_variables); - // Returns true if there is no constraints in this class. - bool IsEmpty() const final { - return num_implications_ == 0 && at_most_ones_.empty(); - } + // Returns the "size" of this class, that is 2 * num_boolean_variables as + // updated by the larger Resize() call. + int64_t literal_size() const { return implications_.size(); } + + // Returns true if no constraints where ever added to this class. + bool IsEmpty() const final { return no_constraint_ever_added_; } // Adds the binary clause (a OR b), which is the same as (not a => b). // Note that it is also equivalent to (not b => a). @@ -721,9 +723,12 @@ class BinaryImplicationGraph : public SatPropagator { // Returns the number of current implications. Note that a => b and not(b) // => not(a) are counted separately since they appear separately in our // propagation lists. The number of size 2 clauses that represent the same - // thing is half this number. - int64_t num_implications() const { return num_implications_; } - int64_t literal_size() const { return implications_.size(); } + // thing is half this number. This should only be used in logs. + int64_t ComputeNumImplicationsForLog() const { + int64_t result = 0; + for (const auto& list : implications_) result += list.size(); + return result; + } // Extract all the binary clauses managed by this class. The Output type must // support an AddBinaryClause(Literal a, Literal b) function. @@ -863,7 +868,7 @@ class BinaryImplicationGraph : public SatPropagator { // // If the final AMO size is smaller than the at_most_one_expansion_size // parameters, we fully expand it. - bool CleanUpAndAddAtMostOnes(int base_index); + ABSL_MUST_USE_RESULT bool CleanUpAndAddAtMostOnes(int base_index); // To be used in DCHECKs(). bool InvariantsAreOk(); @@ -878,6 +883,11 @@ class BinaryImplicationGraph : public SatPropagator { Trail* trail_; DratProofHandler* drat_proof_handler_ = nullptr; + // When problems do not have any implications or at_most_ones this allows to + // reduce the number of work we do here. This will be set to true the first + // time something is added. + bool no_constraint_ever_added_ = true; + // Binary reasons by trail_index. We need a deque because we kept pointers to // elements of this array and this can dynamically change size. std::deque reasons_; diff --git a/ortools/sat/clause_test.cc b/ortools/sat/clause_test.cc index 4a80782cdf..ab70dc266a 100644 --- a/ortools/sat/clause_test.cc +++ b/ortools/sat/clause_test.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include "absl/container/flat_hash_set.h" @@ -194,9 +195,9 @@ TEST(BinaryImplicationGraphTest, TransitiveReduction) { graph->AddBinaryClause(Literal(i, false), Literal(i, true)); } - EXPECT_EQ(graph->num_implications(), 10 * 9); + EXPECT_EQ(graph->ComputeNumImplicationsForLog(), 10 * 9); EXPECT_TRUE(graph->ComputeTransitiveReduction()); - EXPECT_EQ(graph->num_implications(), 9 * 2); + EXPECT_EQ(graph->ComputeNumImplicationsForLog(), 9 * 2); } // This basically just test our DCHECKs. @@ -217,9 +218,9 @@ TEST(BinaryImplicationGraphTest, BasicRandomTransitiveReduction) { graph->AddImplication(Literal(a, true), Literal(b, true)); } - EXPECT_EQ(graph->num_implications(), 2 * num_added); + EXPECT_EQ(graph->ComputeNumImplicationsForLog(), 2 * num_added); EXPECT_TRUE(graph->ComputeTransitiveReduction()); - EXPECT_LT(graph->num_implications(), num_added); + EXPECT_LT(graph->ComputeNumImplicationsForLog(), num_added); } // We generate a random 2-SAT problem, and check that the propagation is @@ -421,6 +422,38 @@ TEST(BinaryImplicationGraphTest, RandomImpliedLiteral) { EXPECT_THAT(seen, UnorderedLiteralsAre(-2, -4, -5, +6, +7)); } +TEST(BinaryImplicationGraphTest, DetectEquivalencePropagateThings) { + Model model; + auto* trail = model.GetOrCreate(); + auto* graph = model.GetOrCreate(); + trail->Resize(10); + graph->Resize(10); + + EXPECT_TRUE(graph->AddAtMostOne(Literals({+1, +2, +3, +4}))); + EXPECT_TRUE(graph->AddAtMostOne(Literals({-2, -1, +3, +4}))); + EXPECT_TRUE(graph->AddAtMostOne(Literals({-4, -1, +2, +3}))); + EXPECT_TRUE(graph->AddAtMostOne(Literals({-3, -1, +2, +4}))); + EXPECT_TRUE(graph->DetectEquivalences()); +} + +void TryAmoEquivalences(absl::Span> cliques) { + Model model; + auto* trail = model.GetOrCreate(); + auto* graph = model.GetOrCreate(); + trail->Resize(1000); + graph->Resize(1000); + for (const auto& clique : cliques) { + std::vector literals; + for (const int i : clique) { + literals.push_back(Literal(i)); + } + if (!graph->AddAtMostOne(literals)) { + return; + } + } + graph->DetectEquivalences(); +} + } // namespace } // namespace sat } // namespace operations_research diff --git a/ortools/sat/constraint_violation.cc b/ortools/sat/constraint_violation.cc index 561173715c..7719cf43a1 100644 --- a/ortools/sat/constraint_violation.cc +++ b/ortools/sat/constraint_violation.cc @@ -53,6 +53,12 @@ int64_t ExprValue(const LinearExpressionProto& expr, return result; } +int64_t AffineValue(const ViewOfAffineLinearExpressionProto& affine, + absl::Span solution) { + if (affine.coeff == 0) return affine.offset; + return affine.coeff * solution[affine.var] + affine.offset; +} + LinearExpressionProto ExprDiff(const LinearExpressionProto& a, const LinearExpressionProto& b) { LinearExpressionProto result; @@ -1238,6 +1244,56 @@ int64_t CompiledNoOverlap2dConstraint::ComputeViolation( return violation; } +template +int64_t CompiledNoOverlap2dWithTwoBoxes::ViolationDelta( + int /*var*/, int64_t /*old_value*/, absl::Span solution) { + if (has_enforcement) { + for (const int lit : enforcements_) { + if (!LiteralValue(lit, solution)) return -violation_; + } + } + + const int64_t x1 = AffineValue(box1_.x_min, solution); + const int64_t X1 = AffineValue(box1_.x_max, solution); + const int64_t x2 = AffineValue(box2_.x_min, solution); + const int64_t X2 = AffineValue(box2_.x_max, solution); + const int64_t repair_x = std::min(X2 - x1, X1 - x2); + if (repair_x <= 0) return -violation_; // disjoint + + const int64_t y1 = AffineValue(box1_.y_min, solution); + const int64_t Y1 = AffineValue(box1_.y_max, solution); + const int64_t y2 = AffineValue(box2_.y_min, solution); + const int64_t Y2 = AffineValue(box2_.y_max, solution); + const int64_t repair_y = std::min(Y2 - y1, Y1 - y2); + if (repair_y <= 0) return -violation_; // disjoint + + const int64_t overlap_x = + std::min(std::max(std::min(X2 - x2, X1 - x1), int64_t{1}), repair_x); + const int64_t overlap_y = + std::min(std::max(std::min(Y2 - y2, Y1 - y1), int64_t{1}), repair_y); + return std::min(repair_x * overlap_y, repair_y * overlap_x) - violation_; +} + +template +std::vector +CompiledNoOverlap2dWithTwoBoxes::UsedVariables( + const CpModelProto& /*model_proto*/) const { + std::vector result; + if (has_enforcement) { + for (const int ref : enforcements_) result.push_back(PositiveRef(ref)); + } + box1_.x_min.AppendVarTo(result); + box1_.x_max.AppendVarTo(result); + box1_.y_min.AppendVarTo(result); + box1_.y_max.AppendVarTo(result); + box2_.x_min.AppendVarTo(result); + box2_.x_max.AppendVarTo(result); + box2_.y_min.AppendVarTo(result); + box2_.y_max.AppendVarTo(result); + gtl::STLSortAndRemoveDuplicates(&result); + return result; +} + // ----- CompiledCircuitConstraint ----- // The violation of a circuit has three parts: @@ -1463,6 +1519,7 @@ LsEvaluator::LsEvaluator(const CpModelProto& cp_model, const SatParameters& params, TimeLimit* time_limit) : cp_model_(cp_model), params_(params), time_limit_(time_limit) { var_to_constraints_.resize(cp_model_.variables_size()); + var_to_dtime_estimate_.resize(cp_model_.variables_size()); jump_value_optimal_.resize(cp_model_.variables_size(), true); num_violated_constraint_per_var_ignoring_objective_.assign( cp_model_.variables_size(), 0); @@ -1481,6 +1538,7 @@ LsEvaluator::LsEvaluator( TimeLimit* time_limit) : cp_model_(cp_model), params_(params), time_limit_(time_limit) { var_to_constraints_.resize(cp_model_.variables_size()); + var_to_dtime_estimate_.resize(cp_model_.variables_size()); jump_value_optimal_.resize(cp_model_.variables_size(), true); num_violated_constraint_per_var_ignoring_objective_.assign( cp_model_.variables_size(), 0); @@ -1498,8 +1556,11 @@ void LsEvaluator::BuildVarConstraintGraph() { for (int ct_index = 0; ct_index < constraints_.size(); ++ct_index) { constraint_to_vars_[ct_index] = constraints_[ct_index]->UsedVariables(cp_model_); + + const double dtime = 1e-8 * constraint_to_vars_[ct_index].size(); for (const int var : constraint_to_vars_[ct_index]) { var_to_constraints_[var].push_back(ct_index); + var_to_dtime_estimate_[var] += dtime; } } @@ -1734,22 +1795,26 @@ void LsEvaluator::CompileOneConstraint(const ConstraintProto& ct) { } for (int i = 0; i + 1 < size; ++i) { - const IntervalConstraintProto& x_interval_i = - cp_model_.constraints(x_intervals[i]).interval(); + const ConstraintProto& x_proto_i = + cp_model_.constraints(x_intervals[i]); + const IntervalConstraintProto& x_interval_i = x_proto_i.interval(); const int64_t x_min_start_i = ExprMin(x_interval_i.start(), cp_model_); const int64_t x_max_end_i = ExprMax(x_interval_i.end(), cp_model_); - const IntervalConstraintProto& y_interval_i = - cp_model_.constraints(y_intervals[i]).interval(); + const ConstraintProto& y_proto_i = + cp_model_.constraints(y_intervals[i]); + const IntervalConstraintProto& y_interval_i = y_proto_i.interval(); const int64_t y_min_start_i = ExprMin(y_interval_i.start(), cp_model_); const int64_t y_max_end_i = ExprMax(y_interval_i.end(), cp_model_); for (int j = i + 1; j < size; ++j) { - const IntervalConstraintProto& x_interval_j = - cp_model_.constraints(x_intervals[j]).interval(); + const ConstraintProto& x_proto_j = + cp_model_.constraints(x_intervals[j]); + const IntervalConstraintProto& x_interval_j = x_proto_j.interval(); const int64_t x_min_start_j = ExprMin(x_interval_j.start(), cp_model_); const int64_t x_max_end_j = ExprMax(x_interval_j.end(), cp_model_); - const IntervalConstraintProto& y_interval_j = - cp_model_.constraints(y_intervals[j]).interval(); + const ConstraintProto& y_proto_j = + cp_model_.constraints(y_intervals[j]); + const IntervalConstraintProto& y_interval_j = y_proto_j.interval(); const int64_t y_min_start_j = ExprMin(y_interval_j.start(), cp_model_); const int64_t y_max_end_j = ExprMax(y_interval_j.end(), cp_model_); @@ -1757,14 +1822,20 @@ void LsEvaluator::CompileOneConstraint(const ConstraintProto& ct) { y_min_start_i >= y_max_end_j || y_min_start_j >= y_max_end_i) { continue; } - ConstraintProto* diffn = expanded_constraints_.add_constraints(); - diffn->mutable_no_overlap_2d()->add_x_intervals(x_intervals[i]); - diffn->mutable_no_overlap_2d()->add_x_intervals(x_intervals[j]); - diffn->mutable_no_overlap_2d()->add_y_intervals(y_intervals[i]); - diffn->mutable_no_overlap_2d()->add_y_intervals(y_intervals[j]); - CompiledNoOverlap2dConstraint* no_overlap_2d = - new CompiledNoOverlap2dConstraint(*diffn, cp_model_); - constraints_.emplace_back(no_overlap_2d); + + const bool has_enforcement = + !x_proto_i.enforcement_literal().empty() || + !x_proto_j.enforcement_literal().empty() || + !y_proto_i.enforcement_literal().empty() || + !y_proto_j.enforcement_literal().empty(); + if (has_enforcement) { + constraints_.emplace_back(new CompiledNoOverlap2dWithTwoBoxes( + x_proto_i, y_proto_i, x_proto_j, y_proto_j)); + } else { + constraints_.emplace_back( + new CompiledNoOverlap2dWithTwoBoxes( + x_proto_i, y_proto_i, x_proto_j, y_proto_j)); + } } } break; @@ -1979,17 +2050,18 @@ double LsEvaluator::WeightedViolationDelta( const int64_t old_value = mutable_solution[var]; mutable_solution[var] += delta; - const int num_linear_constraints = linear_evaluator_.num_constraints(); - for (const int ct_index : var_to_constraints_[var]) { - // We assume linear time delta computation in number of variables. - // TODO(user): refine on a per constraint basis. - dtime_ += 1e-8 * static_cast(constraint_to_vars_[ct_index].size()); + // We assume linear time delta computation in number of variables. + // TODO(user): refine on a per constraint basis. + dtime_ += var_to_dtime_estimate_[var]; + const int num_linear_constraints = linear_evaluator_.num_constraints(); + const std::unique_ptr* data = constraints_.data(); + const auto non_linear_weights = weights.subspan(num_linear_constraints); + for (const int ct_index : var_to_constraints_[var]) { DCHECK_LT(ct_index, constraints_.size()); - const int64_t ct_delta = constraints_[ct_index]->ViolationDelta( - var, old_value, mutable_solution); - result += static_cast(ct_delta) * - weights[ct_index + num_linear_constraints]; + const int64_t ct_delta = + data[ct_index]->ViolationDelta(var, old_value, mutable_solution); + result += static_cast(ct_delta) * non_linear_weights[ct_index]; } // Restore. diff --git a/ortools/sat/constraint_violation.h b/ortools/sat/constraint_violation.h index 014ffab1fc..a77ca1b6a0 100644 --- a/ortools/sat/constraint_violation.h +++ b/ortools/sat/constraint_violation.h @@ -23,6 +23,7 @@ #include "absl/container/flat_hash_map.h" #include "absl/types/span.h" +#include "ortools/base/stl_util.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/util.h" @@ -444,6 +445,7 @@ class LsEvaluator { LinearIncrementalEvaluator linear_evaluator_; std::vector> constraints_; std::vector> var_to_constraints_; + std::vector var_to_dtime_estimate_; std::vector> constraint_to_vars_; std::vector jump_value_optimal_; TimeLimit* time_limit_; @@ -572,8 +574,8 @@ class NoOverlapBetweenTwoIntervals : public CompiledConstraint { class CompiledNoOverlap2dConstraint : public CompiledConstraintWithProto { public: - explicit CompiledNoOverlap2dConstraint(const ConstraintProto& ct_proto, - const CpModelProto& cp_model); + CompiledNoOverlap2dConstraint(const ConstraintProto& ct_proto, + const CpModelProto& cp_model); ~CompiledNoOverlap2dConstraint() override = default; int64_t ComputeViolation(absl::Span solution) override; @@ -582,6 +584,89 @@ class CompiledNoOverlap2dConstraint : public CompiledConstraintWithProto { const CpModelProto& cp_model_; }; +// This is more compact and faster to destroy than storing a +// LinearExpressionProto. +struct ViewOfAffineLinearExpressionProto { + explicit ViewOfAffineLinearExpressionProto( + const LinearExpressionProto& proto) { + if (!proto.vars().empty()) { + DCHECK_EQ(proto.vars().size(), 1); + var = proto.vars(0); + coeff = proto.coeffs(0); + } + offset = proto.offset(); + } + + void AppendVarTo(std::vector& result) const { + if (coeff != 0) result.push_back(var); + } + + int var = 0; + int64_t coeff = 0; + int64_t offset = 0; +}; + +template +class CompiledNoOverlap2dWithTwoBoxes : public CompiledConstraint { + public: + struct Box { + Box(const ConstraintProto& x, const ConstraintProto& y) + : x_min(x.interval().start()), + x_max(x.interval().end()), + y_min(y.interval().start()), + y_max(y.interval().end()) {} + ViewOfAffineLinearExpressionProto x_min; + ViewOfAffineLinearExpressionProto x_max; + ViewOfAffineLinearExpressionProto y_min; + ViewOfAffineLinearExpressionProto y_max; + }; + + CompiledNoOverlap2dWithTwoBoxes(const ConstraintProto& x1, + const ConstraintProto& y1, + const ConstraintProto& x2, + const ConstraintProto& y2) + : box1_(x1, y1), box2_(x2, y2) { + if (has_enforcement) { + enforcements_.assign(x1.enforcement_literal().begin(), + x1.enforcement_literal().end()); + enforcements_.insert(enforcements_.end(), + y1.enforcement_literal().begin(), + y1.enforcement_literal().end()); + enforcements_.insert(enforcements_.end(), + x2.enforcement_literal().begin(), + x2.enforcement_literal().end()); + enforcements_.insert(enforcements_.end(), + y2.enforcement_literal().begin(), + y2.enforcement_literal().end()); + gtl::STLSortAndRemoveDuplicates(&enforcements_); + } + } + + ~CompiledNoOverlap2dWithTwoBoxes() final = default; + + int64_t ComputeViolation(absl::Span solution) final { + // Optimization hack: If we create a ComputeViolationInternal() that we call + // from here and in ViolationDelta(), then the later is not inlined below in + // ViolationDelta() where it matter a lot for performance. + violation_ = 0; + violation_ = ViolationDelta(0, 0, solution); + return violation_; + } + + // Note(user): this is the same implementation as the base one, but it + // avoid one virtual call ! + int64_t ViolationDelta( + int /*var*/, int64_t /*old_value*/, + absl::Span solution_with_new_value) final; + + std::vector UsedVariables(const CpModelProto& model_proto) const final; + + private: + std::vector enforcements_; + const Box box1_; + const Box box2_; +}; + // This can be used to encode reservoir or a cumulative constraints for LS. We // have a set of event time, and we use for overal violation the sum of overload // over time. diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index 6b2b6568cc..c7aa86ffba 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -638,7 +638,8 @@ bool NonOverlappingRectanglesDisjunctivePropagator:: // Optimization: we only initialize the set if we don't have all task here. absl::flat_hash_set requested_boxes_set; - if (requested_boxes.size() != helper_->NumBoxes()) { + const bool not_all_boxes = requested_boxes.size() != helper_->NumBoxes(); + if (not_all_boxes) { requested_boxes_set = {requested_boxes.begin(), requested_boxes.end()}; } @@ -648,10 +649,7 @@ bool NonOverlappingRectanglesDisjunctivePropagator:: auto fixed_boxes = already_checked_fixed_boxes_.view(); for (int i = temp.size(); --i >= 0;) { const int box = temp[i].task_index; - if (requested_boxes.size() != helper_->NumBoxes() && - !requested_boxes_set.contains(box)) { - continue; - } + if (not_all_boxes && !requested_boxes_set.contains(box)) continue; // By definition, fixed boxes are always present. // Doing this check optimize a bit the case where we have many fixed boxes. @@ -910,43 +908,36 @@ bool RectanglePairwisePropagator::Propagate() { point_zero_area_boxes_.clear(); fixed_non_zero_area_boxes_.clear(); non_fixed_non_zero_area_boxes_.clear(); - fixed_non_zero_area_rectangles_.clear(); for (int b : helper_->connected_components()[component_index]) { if (!helper_->IsPresent(b)) continue; const auto [x_size_max, y_size_max] = helper_->GetBoxSizesMax(b); - ItemWithVariableSize* box; + ItemWithVariableSize box = helper_->GetItemWithVariableSize(b); if (x_size_max == 0) { if (y_size_max == 0) { - box = &point_zero_area_boxes_.emplace_back(); + point_zero_area_boxes_.push_back(std::move(box)); } else { - box = &vertical_zero_area_boxes_.emplace_back(); + vertical_zero_area_boxes_.push_back(std::move(box)); } } else if (y_size_max == 0) { - box = &horizontal_zero_area_boxes_.emplace_back(); + horizontal_zero_area_boxes_.push_back(std::move(box)); } else { - if (helper_->IsFixed(b)) { - box = &fixed_non_zero_area_boxes_.emplace_back(); - fixed_non_zero_area_rectangles_.push_back( - helper_->GetItemRangeForSizeMin(b).bounding_area); + // TODO(user): if the number of fixed boxes is large, we keep + // reconstructing them each time this is called for no reason. + if (box.IsFixed()) { + fixed_non_zero_area_boxes_.push_back(std::move(box)); } else { - box = &non_fixed_non_zero_area_boxes_.emplace_back(); + non_fixed_non_zero_area_boxes_.push_back(std::move(box)); } } - *box = helper_->GetItemWithVariableSize(b); } // We ignore pairs of two fixed boxes. The only thing to propagate between // two fixed boxes is a conflict and it should already have been taken care // of by the MandatoryOverlapPropagator propagator. - RETURN_IF_FALSE(FindRestrictionsAndPropagateConflict( non_fixed_non_zero_area_boxes_, fixed_non_zero_area_boxes_, &restrictions)); - RETURN_IF_FALSE(FindRestrictionsAndPropagateConflict( - fixed_non_zero_area_boxes_, non_fixed_non_zero_area_boxes_, - &restrictions)); - // Check zero area boxes against non-zero area boxes. for (auto& non_zero_area_boxes : {fixed_non_zero_area_boxes_, non_fixed_non_zero_area_boxes_}) { @@ -968,24 +959,6 @@ bool RectanglePairwisePropagator::Propagate() { return true; } -bool RectanglePairwisePropagator::FindRestrictionsAndPropagateConflict( - absl::Span items, - std::vector* restrictions) { - const int max_pairs = - params_->max_pairs_pairwise_reasoning_in_no_overlap_2d(); - if (items.size() * (items.size() - 1) / 2 > max_pairs) { - return true; - } - AppendPairwiseRestrictions(items, restrictions); - for (const PairwiseRestriction& restriction : *restrictions) { - if (restriction.type == - PairwiseRestriction::PairwiseRestrictionType::CONFLICT) { - RETURN_IF_FALSE(PropagateTwoBoxes(restriction)); - } - } - return true; -} - bool RectanglePairwisePropagator::FindRestrictionsAndPropagateConflict( absl::Span items1, absl::Span items2, diff --git a/ortools/sat/diffn.h b/ortools/sat/diffn.h index e517473b41..b0fe1b09cb 100644 --- a/ortools/sat/diffn.h +++ b/ortools/sat/diffn.h @@ -167,10 +167,6 @@ class RectanglePairwisePropagator : public PropagatorInterface { delete; // Return false if a conflict is found. - bool FindRestrictionsAndPropagateConflict( - absl::Span items, - std::vector* restrictions); - bool FindRestrictionsAndPropagateConflict( absl::Span items1, absl::Span items2, @@ -186,7 +182,6 @@ class RectanglePairwisePropagator : public PropagatorInterface { int64_t num_pairwise_conflicts_ = 0; int64_t num_pairwise_propagations_ = 0; - std::vector fixed_non_zero_area_rectangles_; std::vector fixed_non_zero_area_boxes_; std::vector non_fixed_non_zero_area_boxes_; std::vector horizontal_zero_area_boxes_; diff --git a/ortools/sat/diffn_util.cc b/ortools/sat/diffn_util.cc index f56ada77c6..ad2f054b93 100644 --- a/ortools/sat/diffn_util.cc +++ b/ortools/sat/diffn_util.cc @@ -599,55 +599,51 @@ void AppendPairwiseRestriction(const ItemWithVariableSize& item1, return; } + PairwiseRestriction::PairwiseRestrictionType type; switch (state) { case 0: // Conflict. The two boxes must overlap in both dimensions. - result->push_back( - {.first_index = item1.index, - .second_index = item2.index, - .type = PairwiseRestriction::PairwiseRestrictionType::CONFLICT}); + type = PairwiseRestriction::PairwiseRestrictionType::CONFLICT; break; case 1: // box2 can only be after box1 on x. if (item1.x.end_min > item2.x.start_min || item2.x.start_max < item1.x.end_max) { - result->push_back({.first_index = item1.index, - .second_index = item2.index, - .type = PairwiseRestriction:: - PairwiseRestrictionType::FIRST_LEFT_OF_SECOND}); + type = + PairwiseRestriction::PairwiseRestrictionType::FIRST_LEFT_OF_SECOND; + break; } - break; - case 2: // box1 an only be after box2 on x. + return; + case 2: + // box1 can only be after box2 on x. if (item2.x.end_min > item1.x.start_min || item1.x.start_max < item2.x.end_max) { - result->push_back({.first_index = item1.index, - .second_index = item2.index, - .type = PairwiseRestriction:: - PairwiseRestrictionType::FIRST_RIGHT_OF_SECOND}); + type = + PairwiseRestriction::PairwiseRestrictionType::FIRST_RIGHT_OF_SECOND; + break; } - break; + return; case 4: - // box2 an only be after box1 on y. + // box2 can only be after box1 on y. if (item1.y.end_min > item2.y.start_min || item2.y.start_max < item1.y.end_max) { - result->push_back({.first_index = item1.index, - .second_index = item2.index, - .type = PairwiseRestriction:: - PairwiseRestrictionType::FIRST_BELOW_SECOND}); + type = PairwiseRestriction::PairwiseRestrictionType::FIRST_BELOW_SECOND; + break; } - break; - case 8: // box1 an only be after box2 on y. + return; + case 8: + // box1 can only be after box2 on y. if (item2.y.end_min > item1.y.start_min || item1.y.start_max < item2.y.end_max) { - result->push_back({.first_index = item1.index, - .second_index = item2.index, - .type = PairwiseRestriction:: - PairwiseRestrictionType::FIRST_ABOVE_SECOND}); + type = PairwiseRestriction::PairwiseRestrictionType::FIRST_ABOVE_SECOND; + break; } - break; + return; default: - break; + return; } + result->push_back( + {.first_index = item1.index, .second_index = item2.index, .type = type}); } } // namespace @@ -1942,8 +1938,7 @@ PostProcessedResult ConvertToRectangle32WithNonZeroSizes( IntegerValue prev_y = 0; Event prev_event = Event::kEnd; int cur_index = -1; - for (int i = 0; i < y_events.size(); ++i) { - const auto [y, event, index] = y_events[i]; + for (const auto [y, event, index] : y_events) { if ((prev_event != event && prev_event != Event::kEnd) || prev_y != y || event == Event::kPoint || cur_index == -1) { ++cur_index; @@ -1974,8 +1969,7 @@ PostProcessedResult ConvertToRectangle32WithNonZeroSizes( IntegerValue prev_x = 0; prev_event = Event::kEnd; cur_index = -1; - for (int i = 0; i < x_events.size(); ++i) { - const auto [x, event, index] = x_events[i]; + for (const auto [x, event, index] : x_events) { if ((prev_event != event && prev_event != Event::kEnd) || prev_x != x || event == Event::kPoint || cur_index == -1) { ++cur_index; @@ -2000,8 +1994,7 @@ PostProcessedResult ConvertToRectangle32WithNonZeroSizes( std::vector sorted_rectangles32; sorted_rectangles32.reserve(rectangles.size()); - for (int i = 0; i < x_events.size(); ++i) { - const auto [x, event, index] = x_events[i]; + for (const auto [x, event, index] : x_events) { if (event == Event::kBegin || event == Event::kPoint) { sorted_rectangles32.push_back(rectangles32[index]); } diff --git a/ortools/sat/diffn_util.h b/ortools/sat/diffn_util.h index 601ce37a73..5cef31ac91 100644 --- a/ortools/sat/diffn_util.h +++ b/ortools/sat/diffn_util.h @@ -269,15 +269,18 @@ std::vector GetIntervalArticulationPoints( std::vector* intervals); struct ItemWithVariableSize { - int index; struct Interval { + bool IsFixed() const { + return start_min == start_max && end_min == end_max; + } + IntegerValue start_min; IntegerValue start_max; IntegerValue end_min; IntegerValue end_max; }; - Interval x; - Interval y; + + bool IsFixed() const { return x.IsFixed() && y.IsFixed(); } template friend void AbslStringify(Sink& sink, const ItemWithVariableSize& item) { @@ -286,6 +289,10 @@ struct ItemWithVariableSize { item.x.end_max, item.y.start_min, item.y.start_max, item.y.end_min, item.y.end_max); } + + int index; + Interval x; + Interval y; }; struct PairwiseRestriction { diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index c217d0fb8c..c74909f5fd 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -1524,6 +1524,7 @@ bool FixedModuloPropagator::PropagateSignsAndTargetRange() { // The sign of target_ is fixed by the sign of expr_. if (integer_trail_->LowerBound(expr_) >= 0 && integer_trail_->LowerBound(target_) < 0) { + // expr >= 0 => target >= 0. if (!integer_trail_->SafeEnqueue(target_.GreaterOrEqual(0), {expr_.GreaterOrEqual(0)})) { return false; @@ -1532,12 +1533,31 @@ bool FixedModuloPropagator::PropagateSignsAndTargetRange() { if (integer_trail_->UpperBound(expr_) <= 0 && integer_trail_->UpperBound(target_) > 0) { + // expr <= 0 => target <= 0. if (!integer_trail_->SafeEnqueue(target_.LowerOrEqual(0), {expr_.LowerOrEqual(0)})) { return false; } } + if (integer_trail_->LowerBound(target_) > 0 && + integer_trail_->LowerBound(expr_) <= 0) { + // target > 0 => expr > 0. + if (!integer_trail_->SafeEnqueue(expr_.GreaterOrEqual(1), + {target_.GreaterOrEqual(1)})) { + return false; + } + } + + if (integer_trail_->UpperBound(target_) < 0 && + integer_trail_->UpperBound(expr_) >= 0) { + // target < 0 => expr < 0. + if (!integer_trail_->SafeEnqueue(expr_.LowerOrEqual(-1), + {target_.LowerOrEqual(-1)})) { + return false; + } + } + return true; } diff --git a/ortools/sat/opb_reader.h b/ortools/sat/opb_reader.h index a6fb2a4441..8be4d750e1 100644 --- a/ortools/sat/opb_reader.h +++ b/ortools/sat/opb_reader.h @@ -16,14 +16,22 @@ #include #include +#include #include +#include #include +#include "absl/container/flat_hash_map.h" +#include "absl/container/flat_hash_set.h" #include "absl/log/check.h" +#include "absl/log/log.h" #include "absl/strings/numbers.h" #include "absl/strings/str_split.h" +#include "absl/types/span.h" #include "ortools/base/logging.h" -#include "ortools/sat/boolean_problem.pb.h" +#include "ortools/base/stl_util.h" +#include "ortools/sat/cp_model.pb.h" +#include "ortools/sat/cp_model_utils.h" #include "ortools/util/filelineiter.h" namespace operations_research { @@ -31,7 +39,7 @@ namespace sat { // This class loads a file in pbo file format into a LinearBooleanProblem. // The format is described here: -// http://www.cril.univ-artois.fr/PB12/format.pdf +// http://www.cril.univ-artois.fr/PB24/format.pdf class OpbReader { public: OpbReader() = default; @@ -39,25 +47,69 @@ class OpbReader { OpbReader(const OpbReader&) = delete; OpbReader& operator=(const OpbReader&) = delete; + // Returns the number of variables in the problem. + int num_variables() const { return num_variables_; } + // Loads the given opb filename into the given problem. - bool Load(const std::string& filename, LinearBooleanProblem* problem) { - problem->Clear(); - problem->set_name(ExtractProblemName(filename)); + // Returns true on success. + bool LoadAndValidate(const std::string& filename, CpModelProto* model) { + model->Clear(); + model->set_name(ExtractProblemName(filename)); num_variables_ = 0; int num_lines = 0; + + // Read constraints line by line (1 constraint per line). + // We process into a temporary structure to support non linear constraints + // and weighted constraints. for (const std::string& line : FileLines(filename)) { ++num_lines; - ProcessNewLine(problem, line); + ProcessNewLine(line); } if (num_lines == 0) { - LOG(FATAL) << "File '" << filename << "' is empty or can't be read."; + LOG(ERROR) << "File '" << filename << "' is empty or can't be read."; + return false; } - problem->set_num_variables(num_variables_); + + LOG(INFO) << "Read " << num_lines << " lines from " << filename; + LOG(INFO) << "#variables: " << num_variables_; + LOG(INFO) << "#constraints: " << constraints_.size(); + LOG(INFO) << "#objective: " << objective_.size(); + + const std::string error_message = ValidateModel(); + if (!error_message.empty()) { + LOG(ERROR) << "Error while trying to parse '" << filename + << "': " << error_message; + return false; + } + + BuildModel(model); return true; } private: + // A term is coeff * Product(literals). + // Note that it is okay to have duplicate literals here, we will just merge + // them. Having a literal and its negation will always result in a product of + // zero. + struct PbTerm { + int64_t coeff; + std::vector literals; // CpModelProto literals + }; + + enum PbConstraintType { + UNDEFINED_OPERATION, + GE_OPERATION, + EQ_OPERATION, + }; + + struct PbConstraint { + std::vector terms; + PbConstraintType type = UNDEFINED_OPERATION; + int64_t rhs = std::numeric_limits::min(); + int64_t soft_cost = std::numeric_limits::max(); + }; + // Since the problem name is not stored in the cnf format, we infer it from // the file name. static std::string ExtractProblemName(const std::string& filename) { @@ -67,71 +119,233 @@ class OpbReader { return problem_name; } - void ProcessNewLine(LinearBooleanProblem* problem, const std::string& line) { + void ProcessNewLine(const std::string& line) { const std::vector words = absl::StrSplit(line, absl::ByAnyChar(" ;"), absl::SkipEmpty()); if (words.empty() || words[0].empty() || words[0][0] == '*') { + // TODO(user): Parse comments. return; } if (words[0] == "min:") { - LinearObjective* objective = problem->mutable_objective(); for (int i = 1; i < words.size(); ++i) { const std::string& word = words[i]; if (word.empty() || word[0] == ';') continue; if (word[0] == 'x') { - int literal; - CHECK(absl::SimpleAtoi(word.substr(1), &literal)); - num_variables_ = std::max(num_variables_, literal); - objective->add_literals(literal); + int index; + CHECK(absl::SimpleAtoi(word.substr(1), &index)); + num_variables_ = std::max(num_variables_, index); + objective_.back().literals.push_back( + PbLiteralToCpModelLiteral(index)); + } else if (word[0] == '~' && word[1] == 'x') { + int index; + CHECK(absl::SimpleAtoi(word.substr(2), &index)); + num_variables_ = std::max(num_variables_, index); + objective_.back().literals.push_back( + NegatedRef(PbLiteralToCpModelLiteral(index))); } else { - int64_t value; - CHECK(absl::SimpleAtoi(word, &value)); - objective->add_coefficients(value); + int64_t coefficient; + CHECK(absl::SimpleAtoi(word, &coefficient)); + objective_.push_back({coefficient, {}}); } } - if (objective->literals_size() != objective->coefficients_size()) { - LOG(INFO) << "words.size() = " << words.size(); - LOG(FATAL) << "Failed to parse objective:\n " << line; + + // Normalize objective literals. + for (PbTerm& term : objective_) { + if (term.literals.size() <= 1) continue; + gtl::STLSortAndRemoveDuplicates(&term.literals); + CHECK_GT(term.literals.size(), 1); } + return; } - LinearBooleanConstraint* constraint = problem->add_constraints(); + + PbConstraint constraint; for (int i = 0; i < words.size(); ++i) { const std::string& word = words[i]; CHECK(!word.empty()); - if (word == ">=") { + if (word[0] == '[') { // Soft constraint. + int64_t soft_cost; + CHECK(absl::SimpleAtoi(word.substr(1, word.size() - 2), &soft_cost)); + constraint.soft_cost = soft_cost; + } else if (word == ">=") { CHECK_LT(i + 1, words.size()); - int64_t value; - CHECK(absl::SimpleAtoi(words[i + 1], &value)); - constraint->set_lower_bound(value); + int64_t rhs; + CHECK(absl::SimpleAtoi(words[i + 1], &rhs)); + constraint.type = GE_OPERATION; + constraint.rhs = rhs; break; } else if (word == "=") { CHECK_LT(i + 1, words.size()); - int64_t value; - CHECK(absl::SimpleAtoi(words[i + 1], &value)); - constraint->set_upper_bound(value); - constraint->set_lower_bound(value); + int64_t rhs; + CHECK(absl::SimpleAtoi(words[i + 1], &rhs)); + constraint.type = EQ_OPERATION; + constraint.rhs = rhs; break; + } else if (word[0] == 'x') { + int index; + CHECK(absl::SimpleAtoi(word.substr(1), &index)); + num_variables_ = std::max(num_variables_, index); + constraint.terms.back().literals.push_back( + PbLiteralToCpModelLiteral(index)); + } else if (word[0] == '~' && word[1] == 'x') { + int index; + CHECK(absl::SimpleAtoi(word.substr(2), &index)); + num_variables_ = std::max(num_variables_, index); + constraint.terms.back().literals.push_back( + NegatedRef(PbLiteralToCpModelLiteral(index))); } else { - if (word[0] == 'x') { - int literal; - CHECK(absl::SimpleAtoi(word.substr(1), &literal)); - num_variables_ = std::max(num_variables_, literal); - constraint->add_literals(literal); - } else { - int64_t value; - CHECK(absl::SimpleAtoi(words[i], &value)); - constraint->add_coefficients(value); + // Note that coefficient always appear before the variable/variables. + int64_t coefficient; + CHECK(absl::SimpleAtoi(word, &coefficient)); + constraint.terms.push_back({coefficient, {}}); + } + } + + // Normalize literals. + for (PbTerm& term : constraint.terms) { + if (term.literals.size() <= 1) continue; + gtl::STLSortAndRemoveDuplicates(&term.literals); + CHECK_GT(term.literals.size(), 1); + } + + constraints_.push_back(std::move(constraint)); + } + + std::string ValidateModel() { + // Normalize and validate constraints. + for (const PbConstraint& constraint : constraints_) { + if (constraint.rhs == std::numeric_limits::min()) { + return "constraint error: undefined rhs"; + } + + if (constraint.type == UNDEFINED_OPERATION) { + return "constraint error: undefined operation"; + } + + for (const PbTerm& term : constraint.terms) { + if (term.coeff == 0) { + return "constraint error: coefficient cannot be zero"; + } + if (term.literals.empty()) return "constraint error: empty literals"; + if (term.literals.size() == 1) { + if (!RefIsPositive(term.literals[0])) { + return "constraint error: linear terms must use positive literals"; + } } } } - if (constraint->literals_size() != constraint->coefficients_size()) { - LOG(FATAL) << "Failed to parse constraint:\n " << line; + + // Normalize and validate objective. + if (objective_.empty()) return ""; // No objective. + + for (const PbTerm& term : objective_) { + if (term.coeff == 0) return "objective error: coefficient cannot be zero"; + if (term.literals.empty()) return "objective error: empty literals"; + if (term.literals.size() == 1) { + if (!RefIsPositive(term.literals[0])) { + return "objective error: linear terms must use positive literals"; + } + return ""; + } + } + + return ""; + } + + static int PbLiteralToCpModelLiteral(int pb_literal) { + return pb_literal > 0 ? pb_literal - 1 : -pb_literal; + } + + int GetVariable(const PbTerm& term, CpModelProto* model) { + CHECK(!term.literals.empty()); + if (term.literals.size() == 1) { + CHECK(RefIsPositive(term.literals[0])); + return term.literals[0]; + } + + const auto it = product_to_var_.find(term.literals); + if (it != product_to_var_.end()) { + return it->second; + } + + const int var_index = model->variables_size(); + IntegerVariableProto* var_proto = model->add_variables(); + var_proto->add_domain(0); + var_proto->add_domain(1); + + product_to_var_[term.literals] = var_index; + + // Link the new variable to the terms. + // var_index => and(literals). + ConstraintProto* var_to_literals = model->add_constraints(); + var_to_literals->add_enforcement_literal(var_index); + + // and(literals) => var_index. + ConstraintProto* literals_to_var = model->add_constraints(); + literals_to_var->mutable_bool_and()->add_literals(var_index); + + for (const int proto_literal : term.literals) { + var_to_literals->mutable_bool_and()->add_literals(proto_literal); + literals_to_var->add_enforcement_literal(proto_literal); + } + + return var_index; + } + + void BuildModel(CpModelProto* model) { + // We know how many variables we have, so we can add them all. + for (int i = 0; i < num_variables_; ++i) { + IntegerVariableProto* var = model->add_variables(); + var->add_domain(0); + var->add_domain(1); + } + + for (const PbConstraint& constraint : constraints_) { + ConstraintProto* ct = model->add_constraints(); + LinearConstraintProto* lin = ct->mutable_linear(); + for (const PbTerm& term : constraint.terms) { + lin->add_vars(GetVariable(term, model)); + lin->add_coeffs(term.coeff); + } + if (constraint.type == GE_OPERATION) { + lin->add_domain(constraint.rhs); + lin->add_domain(std::numeric_limits::max()); + } else if (constraint.type == EQ_OPERATION) { + lin->add_domain(constraint.rhs); + lin->add_domain(constraint.rhs); + } else { + LOG(FATAL) << "Unsupported operation: " << constraint.type; + } + + if (constraint.soft_cost != std::numeric_limits::max()) { + const int violation_var_index = model->variables_size(); + IntegerVariableProto* violation_var = model->add_variables(); + violation_var->add_domain(0); + violation_var->add_domain(1); + + // Update the objective. + model->mutable_objective()->add_vars(violation_var_index); + model->mutable_objective()->add_coeffs(constraint.soft_cost); + + // Add the enforcement literal to ct. + ct->add_enforcement_literal(NegatedRef(violation_var_index)); + } + } + + if (!objective_.empty()) { + CpObjectiveProto* obj = model->mutable_objective(); + for (const PbTerm& term : objective_) { + obj->add_vars(GetVariable(term, model)); + obj->add_coeffs(term.coeff); + } } } int num_variables_; + std::vector objective_; + std::vector constraints_; + absl::flat_hash_map, int> product_to_var_; }; } // namespace sat diff --git a/ortools/sat/sat_inprocessing.cc b/ortools/sat/sat_inprocessing.cc index 3c5ab4257e..796e7ae4aa 100644 --- a/ortools/sat/sat_inprocessing.cc +++ b/ortools/sat/sat_inprocessing.cc @@ -163,7 +163,7 @@ bool Inprocessing::PresolveLoop(SatPresolveOptions options) { logger_, "[Pure SAT presolve]", " num_fixed: ", trail_->Index(), " num_redundant: ", implication_graph_->num_redundant_literals() / 2, "/", sat_solver_->NumVariables(), - " num_implications: ", implication_graph_->num_implications(), + " num_implications: ", implication_graph_->ComputeNumImplicationsForLog(), " num_watched_clauses: ", clause_manager_->num_watched_clauses(), " dtime: ", time_limit_->GetElapsedDeterministicTime() - start_dtime, "/", options.deterministic_time_limit, " wtime: ", wall_timer.Get(), @@ -280,7 +280,7 @@ bool Inprocessing::InprocessingRound() { logger_, "Inprocessing.", " fixed:", trail_->Index(), " equiv:", implication_graph_->num_redundant_literals() / 2, " bools:", sat_solver_->NumVariables(), - " implications:", implication_graph_->num_implications(), + " implications:", implication_graph_->ComputeNumImplicationsForLog(), " watched:", clause_manager_->num_watched_clauses(), " minimization:", mini_num_clause, "|", mini_num_removed, " dtime:", time_limit_->GetElapsedDeterministicTime() - start_dtime, @@ -664,8 +664,7 @@ bool StampingSimplifier::DoOneRound(bool log_info) { num_removed_literals_ = 0; num_fixed_ = 0; - if (implication_graph_->literal_size() == 0) return true; - if (implication_graph_->num_implications() == 0) return true; + if (implication_graph_->IsEmpty()) return true; if (!stamps_are_already_computed_) { // We need a DAG so that we don't have cycle while we sample the tree. @@ -696,8 +695,7 @@ bool StampingSimplifier::ComputeStampsForNextRound(bool log_info) { dtime_ = 0.0; num_fixed_ = 0; - if (implication_graph_->literal_size() == 0) return true; - if (implication_graph_->num_implications() == 0) return true; + if (implication_graph_->IsEmpty()) return true; implication_graph_->RemoveFixedVariables(); if (!implication_graph_->DetectEquivalences(log_info)) return true; diff --git a/ortools/sat/sat_inprocessing_test.cc b/ortools/sat/sat_inprocessing_test.cc index 963ac8a7b5..76091f3df3 100644 --- a/ortools/sat/sat_inprocessing_test.cc +++ b/ortools/sat/sat_inprocessing_test.cc @@ -133,7 +133,7 @@ TEST(InprocessingTest, ClauseSubsumptionAndStrengthening) { // We added {+1, +2} and {+1, -3} there. // TODO(user): make sure we don't add twice the implications. auto* implication_graph = model.GetOrCreate(); - EXPECT_EQ(implication_graph->num_implications(), 6); + EXPECT_EQ(implication_graph->ComputeNumImplicationsForLog(), 6); EXPECT_EQ(implication_graph->Implications(Literal(-1)).size(), 3); EXPECT_THAT(implication_graph->Implications(Literal(-1)), ::testing::UnorderedElementsAre(Literal(+2), Literal(+2), @@ -146,7 +146,7 @@ TEST(InprocessingTest, ClauseSubsumptionAndStrengthening) { // Depending on the implication added, we don't get the same clauses. auto* implication_graph = model.GetOrCreate(); - EXPECT_EQ(implication_graph->num_implications(), 2); + EXPECT_EQ(implication_graph->ComputeNumImplicationsForLog(), 2); EXPECT_EQ(implication_graph->Implications(Literal(-1)).size(), 1); if (implication_graph->Implications(Literal(-1))[0] == Literal(+2)) { EXPECT_EQ(all_clauses[2]->AsSpan(), Literals({+2, +6, +1, +1})); diff --git a/ortools/sat/sat_runner.cc b/ortools/sat/sat_runner.cc index 8f21d79c8b..adaec46cef 100644 --- a/ortools/sat/sat_runner.cc +++ b/ortools/sat/sat_runner.cc @@ -29,8 +29,6 @@ #include "ortools/base/helpers.h" #include "ortools/base/options.h" #include "ortools/base/path.h" -#include "ortools/sat/boolean_problem.h" -#include "ortools/sat/boolean_problem.pb.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_solver.h" #include "ortools/sat/cp_model_utils.h" @@ -95,13 +93,12 @@ std::string ExtractName(absl::string_view full_filename) { bool LoadProblem(const std::string& filename, absl::string_view hint_file, absl::string_view domain_file, CpModelProto* cp_model) { if (absl::EndsWith(filename, ".opb") || - absl::EndsWith(filename, ".opb.bz2")) { + absl::EndsWith(filename, ".opb.bz2") || + absl::EndsWith(filename, ".opb.gz")) { OpbReader reader; - LinearBooleanProblem problem; - if (!reader.Load(filename, &problem)) { + if (!reader.LoadAndValidate(filename, cp_model)) { LOG(FATAL) << "Cannot load file '" << filename << "'."; } - *cp_model = BooleanProblemToCpModelproto(problem); } else if (absl::EndsWith(filename, ".cnf") || absl::EndsWith(filename, ".cnf.xz") || absl::EndsWith(filename, ".cnf.gz") || diff --git a/ortools/sat/sat_solver.cc b/ortools/sat/sat_solver.cc index 89a44a9959..69d101f367 100644 --- a/ortools/sat/sat_solver.cc +++ b/ortools/sat/sat_solver.cc @@ -1408,7 +1408,7 @@ SatSolver::Status SatSolver::SolveInternal(TimeLimit* time_limit, SOLVER_LOG(logger_, "Number of clauses (size > 2): ", clauses_propagator_->num_clauses()); SOLVER_LOG(logger_, "Number of binary clauses: ", - binary_implication_graph_->num_implications()); + binary_implication_graph_->ComputeNumImplicationsForLog()); SOLVER_LOG(logger_, "Number of linear constraints: ", pb_constraints_->NumberOfConstraints()); SOLVER_LOG(logger_, "Number of fixed variables: ", trail_->Index()); @@ -1812,7 +1812,8 @@ std::string SatSolver::RunningStatisticsString() const { clauses_propagator_->num_clauses() - clauses_propagator_->num_removable_clauses(), clauses_propagator_->num_removable_clauses(), - binary_implication_graph_->num_implications(), restart_->NumRestarts(), + binary_implication_graph_->ComputeNumImplicationsForLog(), + restart_->NumRestarts(), num_variables_.value() - num_processed_fixed_variables_); } diff --git a/ortools/sat/scheduling_helpers.cc b/ortools/sat/scheduling_helpers.cc index bbca7a12f1..9b10e4c91b 100644 --- a/ortools/sat/scheduling_helpers.cc +++ b/ortools/sat/scheduling_helpers.cc @@ -246,7 +246,9 @@ void SchedulingConstraintHelper::InitSortedVectors() { recompute_all_cache_ = true; recompute_cache_.Resize(num_tasks); + non_fixed_intervals_.resize(num_tasks); for (int t = 0; t < num_tasks; ++t) { + non_fixed_intervals_[t] = t; recompute_cache_.Set(t); } @@ -313,9 +315,23 @@ bool SchedulingConstraintHelper::SynchronizeAndSetTimeDirection( } if (recompute_all_cache_) { - for (int t = 0; t < recompute_cache_.size(); ++t) { + for (const int t : non_fixed_intervals_) { if (!UpdateCachedValues(t)) return false; } + + // We also update non_fixed_intervals_ at level zero so that we will never + // scan them again. + if (sat_solver_->CurrentDecisionLevel() == 0) { + int new_size = 0; + for (const int t : non_fixed_intervals_) { + if (IsPresent(t) && StartIsFixed(t) && EndIsFixed(t) && + SizeIsFixed(t)) { + continue; + } + non_fixed_intervals_[new_size++] = t; + } + non_fixed_intervals_.resize(new_size); + } } else { for (const int t : recompute_cache_) { if (!UpdateCachedValues(t)) return false; diff --git a/ortools/sat/scheduling_helpers.h b/ortools/sat/scheduling_helpers.h index b4353868f8..def922023e 100644 --- a/ortools/sat/scheduling_helpers.h +++ b/ortools/sat/scheduling_helpers.h @@ -457,6 +457,11 @@ class SchedulingConstraintHelper : public PropagatorInterface { bool recompute_all_cache_ = true; Bitset64 recompute_cache_; + // For large problems, LNS will have a lot of fixed intervals. + // And fixed intervals will never changes, so there is no point recomputing + // the cache for them. + std::vector non_fixed_intervals_; + // Reason vectors. std::vector literal_reason_; std::vector integer_reason_;