diff --git a/ortools/sat/2d_distances_propagator.cc b/ortools/sat/2d_distances_propagator.cc new file mode 100644 index 0000000000..62eb603ee8 --- /dev/null +++ b/ortools/sat/2d_distances_propagator.cc @@ -0,0 +1,223 @@ +// Copyright 2010-2025 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/2d_distances_propagator.h" + +#include +#include +#include +#include + +#include "absl/container/flat_hash_map.h" +#include "absl/log/check.h" +#include "absl/log/log.h" +#include "absl/types/span.h" +#include "ortools/base/stl_util.h" +#include "ortools/sat/cp_model.pb.h" +#include "ortools/sat/integer.h" +#include "ortools/sat/integer_base.h" +#include "ortools/sat/linear_propagation.h" +#include "ortools/sat/model.h" +#include "ortools/sat/no_overlap_2d_helper.h" +#include "ortools/sat/precedences.h" +#include "ortools/sat/sat_base.h" +#include "ortools/sat/scheduling_helpers.h" +#include "ortools/sat/synchronization.h" + +namespace operations_research { +namespace sat { + +Precedences2DPropagator::Precedences2DPropagator( + NoOverlap2DConstraintHelper* helper, Model* model) + : helper_(*helper), + binary_relations_maps_(model->GetOrCreate()), + shared_stats_(model->GetOrCreate()) { + model->GetOrCreate()->SetPushAffineUbForBinaryRelation(); +} + +void Precedences2DPropagator::CollectPairsOfBoxesWithNonTrivialDistance() { + helper_.SynchronizeAndSetDirection(); + non_trivial_pairs_.clear(); + + struct VarUsage { + // boxes[0=x, 1=y][0=start, 1=end] + std::vector boxes[2][2]; + }; + absl::flat_hash_map var_to_box_and_coeffs; + + for (int dim = 0; dim < 2; ++dim) { + const SchedulingConstraintHelper& dim_helper = + dim == 0 ? helper_.x_helper() : helper_.y_helper(); + for (int i = 0; i < helper_.NumBoxes(); ++i) { + const absl::Span interval_points = + i == 0 ? dim_helper.Starts() : dim_helper.Ends(); + for (int j = 0; j < 2; ++j) { + if (interval_points[i].var != kNoIntegerVariable) { + var_to_box_and_coeffs[PositiveVariable(interval_points[i].var)] + .boxes[dim][j] + .push_back(i); + } + } + } + } + + VLOG(2) << "CollectPairsOfBoxesWithNonTrivialDistance called, num_exprs: " + << binary_relations_maps_->GetAllExpressionsWithAffineBounds().size(); + for (const LinearExpression2& expr : + binary_relations_maps_->GetAllExpressionsWithAffineBounds()) { + auto it1 = var_to_box_and_coeffs.find(PositiveVariable(expr.vars[0])); + auto it2 = var_to_box_and_coeffs.find(PositiveVariable(expr.vars[1])); + if (it1 == var_to_box_and_coeffs.end() || + it2 == var_to_box_and_coeffs.end()) { + continue; + } + + const VarUsage& usage1 = it1->second; + const VarUsage& usage2 = it2->second; + for (int dim = 0; dim < 2; ++dim) { + const SchedulingConstraintHelper& dim_helper = + dim == 0 ? helper_.x_helper() : helper_.y_helper(); + for (const int box1 : usage1.boxes[dim][0 /* start */]) { + for (const int box2 : usage2.boxes[dim][1 /* end */]) { + const AffineExpression& start = dim_helper.Starts()[box1]; + const AffineExpression& end = dim_helper.Ends()[box2]; + LinearExpression2 expr2; + expr2.vars[0] = start.var; + expr2.vars[1] = NegationOf(end.var); + expr2.coeffs[0] = start.coeff; + expr2.coeffs[1] = end.coeff; + expr2.SimpleCanonicalization(); + expr2.DivideByGcd(); + if (expr == expr2) { + if (box1 < box2) { + non_trivial_pairs_.push_back({box1, box2}); + } else { + non_trivial_pairs_.push_back({box2, box1}); + } + } + } + } + } + } + + gtl::STLSortAndRemoveDuplicates(&non_trivial_pairs_); +} + +bool Precedences2DPropagator::Propagate() { + if (!helper_.SynchronizeAndSetDirection()) return false; + if (last_helper_inprocessing_count_ != helper_.InProcessingCount() || + helper_.x_helper().CurrentDecisionLevel() == 0 || + last_num_expressions_ != + binary_relations_maps_->NumExpressionsWithAffineBounds()) { + last_helper_inprocessing_count_ = helper_.InProcessingCount(); + last_num_expressions_ = + binary_relations_maps_->NumExpressionsWithAffineBounds(); + CollectPairsOfBoxesWithNonTrivialDistance(); + } + + num_calls_++; + + SchedulingConstraintHelper* helpers[2] = {&helper_.x_helper(), + &helper_.y_helper()}; + + for (const auto& [box1, box2] : non_trivial_pairs_) { + DCHECK(box1 < helper_.NumBoxes()); + DCHECK(box2 < helper_.NumBoxes()); + if (!helper_.IsPresent(box1) && !helper_.IsPresent(box2)) { + continue; + } + + bool is_unfeasible = true; + for (int dim = 0; dim < 2; dim++) { + const SchedulingConstraintHelper* helper = helpers[dim]; + for (int j = 0; j < 2; j++) { + int b1 = box1; + int b2 = box2; + if (j == 1) { + std::swap(b1, b2); + } + LinearExpression2 expr; + expr.vars[0] = helper->Starts()[b1].var; + expr.vars[1] = NegationOf(helper->Ends()[b2].var); + expr.coeffs[0] = helper->Starts()[b1].coeff; + expr.coeffs[1] = helper->Ends()[b2].coeff; + const IntegerValue ub_of_start_minus_end_value = + binary_relations_maps_->UpperBound(expr) + + helper->Starts()[b1].constant - helper->Ends()[b2].constant; + if (ub_of_start_minus_end_value >= 0) { + is_unfeasible = false; + break; + } + if (!is_unfeasible) break; + } + } + if (!is_unfeasible) continue; + + // We have a mandatory overlap on both x and y! Explain and propagate. + + helper_.ClearReason(); + num_conflicts_++; + std::vector reason; + std::vector lit_reason; + + for (int dim = 0; dim < 2; dim++) { + SchedulingConstraintHelper* helper = helpers[dim]; + for (int j = 0; j < 2; j++) { + int b1 = box1; + int b2 = box2; + if (j == 1) { + std::swap(b1, b2); + } + LinearExpression2 expr; + expr.vars[0] = helper->Starts()[b1].var; + expr.vars[1] = NegationOf(helper->Ends()[b2].var); + expr.coeffs[0] = helper->Starts()[b1].coeff; + expr.coeffs[1] = helper->Ends()[b2].coeff; + binary_relations_maps_->AddReasonForUpperBoundLowerThan( + expr, + -(helper->Starts()[b1].constant - helper->Ends()[b2].constant), + &lit_reason, &reason); + } + } + helper_.AddPresenceReason(box1); + helper_.AddPresenceReason(box2); + helper_.x_helper().MutableIntegerReason()->insert( + helper_.x_helper().MutableIntegerReason()->end(), reason.begin(), + reason.end()); + helper_.x_helper().MutableLiteralReason()->insert( + helper_.x_helper().MutableLiteralReason()->end(), lit_reason.begin(), + lit_reason.end()); + return helper_.ReportConflict(); + } + return true; +} + +int Precedences2DPropagator::RegisterWith(GenericLiteralWatcher* watcher) { + const int id = watcher->Register(this); + helper_.WatchAllBoxes(id); + // TODO(user): add an API to BinaryRelationsMaps to watch linear2 + return id; +} + +Precedences2DPropagator::~Precedences2DPropagator() { + if (!VLOG_IS_ON(1)) return; + std::vector> stats; + stats.push_back({"Precedences2DPropagator/called", num_calls_}); + stats.push_back({"Precedences2DPropagator/conflicts", num_conflicts_}); + stats.push_back({"Precedences2DPropagator/pairs", non_trivial_pairs_.size()}); + + shared_stats_->AddStats(stats); +} + +} // namespace sat +} // namespace operations_research diff --git a/ortools/sat/2d_distances_propagator.h b/ortools/sat/2d_distances_propagator.h new file mode 100644 index 0000000000..b05e6b1a3d --- /dev/null +++ b/ortools/sat/2d_distances_propagator.h @@ -0,0 +1,67 @@ +// Copyright 2010-2025 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_2D_DISTANCES_PROPAGATOR_H_ +#define OR_TOOLS_SAT_2D_DISTANCES_PROPAGATOR_H_ + +#include +#include +#include + +#include "ortools/sat/integer.h" +#include "ortools/sat/model.h" +#include "ortools/sat/no_overlap_2d_helper.h" +#include "ortools/sat/precedences.h" +#include "ortools/sat/synchronization.h" + +namespace operations_research { +namespace sat { + +// This class implements a propagator for non_overlap_2d constraints that uses +// the BinaryRelationsMaps to detect precedences between pairs of boxes and +// detect a conflict if the precedences implies an overlap between the two +// boxes. For doing this efficiently, it keep track of pairs of boxes that have +// non-fixed precedences in the BinaryRelationsMaps and only check those in the +// propagation. +class Precedences2DPropagator : public PropagatorInterface { + public: + Precedences2DPropagator(NoOverlap2DConstraintHelper* helper, Model* model); + + ~Precedences2DPropagator() override; + + bool Propagate() final; + int RegisterWith(GenericLiteralWatcher* watcher); + + private: + void CollectPairsOfBoxesWithNonTrivialDistance(); + + std::vector> non_trivial_pairs_; + + NoOverlap2DConstraintHelper& helper_; + BinaryRelationsMaps* binary_relations_maps_; + SharedStatistics* shared_stats_; + + int last_helper_inprocessing_count_ = -1; + int last_num_expressions_ = -1; + + int64_t num_conflicts_ = 0; + int64_t num_calls_ = 0; + + Precedences2DPropagator(const Precedences2DPropagator&) = delete; + Precedences2DPropagator& operator=(const Precedences2DPropagator&) = delete; +}; + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_2D_DISTANCES_PROPAGATOR_H_ diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index 33c7851df4..44c42cfcd1 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -114,6 +114,29 @@ proto_library( srcs = ["cp_model.proto"], ) +cc_library( + name = "2d_distances_propagator", + srcs = ["2d_distances_propagator.cc"], + hdrs = ["2d_distances_propagator.h"], + deps = [ + ":cp_model_cc_proto", + ":integer", + ":integer_base", + ":linear_propagation", + ":model", + ":no_overlap_2d_helper", + ":precedences", + ":sat_base", + ":scheduling_helpers", + ":synchronization", + "//ortools/base:stl_util", + "@abseil-cpp//absl/container:flat_hash_map", + "@abseil-cpp//absl/log", + "@abseil-cpp//absl/log:check", + "@abseil-cpp//absl/types:span", + ], +) + cc_library( name = "2d_mandatory_overlap_propagator", srcs = ["2d_mandatory_overlap_propagator.cc"], @@ -809,6 +832,7 @@ cc_library( srcs = ["cp_model_loader.cc"], hdrs = ["cp_model_loader.h"], deps = [ + ":2d_distances_propagator", ":all_different", ":circuit", ":clause", @@ -3175,6 +3199,7 @@ cc_test( "//ortools/util:random_engine", "//ortools/util:sorted_interval_list", "@abseil-cpp//absl/container:btree", + "@abseil-cpp//absl/container:flat_hash_set", "@abseil-cpp//absl/log:check", "@abseil-cpp//absl/numeric:int128", "@abseil-cpp//absl/random", @@ -3514,8 +3539,8 @@ cc_test( "@abseil-cpp//absl/random:distributions", "@abseil-cpp//absl/strings", "@abseil-cpp//absl/types:span", - "@google_benchmark//:benchmark", "@fuzztest//fuzztest:fuzztest_gtest_main", + "@google_benchmark//:benchmark", ], ) @@ -3524,6 +3549,7 @@ cc_library( srcs = ["diffn.cc"], hdrs = ["diffn.h"], deps = [ + ":2d_distances_propagator", ":2d_mandatory_overlap_propagator", ":2d_orthogonal_packing", ":2d_try_edge_propagator", diff --git a/ortools/sat/cp_model.cc b/ortools/sat/cp_model.cc index 098bd85db8..30a888d3d8 100644 --- a/ortools/sat/cp_model.cc +++ b/ortools/sat/cp_model.cc @@ -47,10 +47,10 @@ BoolVar BoolVar::WithName(absl::string_view name) { std::string BoolVar::Name() const { if (builder_ == nullptr) return "null"; - const std::string& name = + absl::string_view name = builder_->Proto().variables(PositiveRef(index_)).name(); if (RefIsPositive(index_)) { - return name; + return std::string(name); } else { return absl::StrCat("Not(", name, ")"); } diff --git a/ortools/sat/cp_model_loader.cc b/ortools/sat/cp_model_loader.cc index 6f5af4cad2..f053299657 100644 --- a/ortools/sat/cp_model_loader.cc +++ b/ortools/sat/cp_model_loader.cc @@ -1271,49 +1271,29 @@ void LoadLinearConstraint(const ConstraintProto& ct, Model* m) { rhs_max = std::min(rhs_max, max_sum.value()); if (vars.size() == 2) { - if (std::abs(coeffs[0]) == std::abs(coeffs[1])) { - const int64_t magnitude = std::abs(coeffs[0]); - IntegerVariable v1 = vars[0]; - IntegerVariable v2 = vars[1]; - if (coeffs[0] < 0) v1 = NegationOf(v1); - if (coeffs[1] > 0) v2 = NegationOf(v2); - - // magnitude * v1 <= magnitude * v2 + rhs_max. - precedences->Add(v1, v2, MathUtil::CeilOfRatio(-rhs_max, magnitude)); - - // magnitude * v1 >= magnitude * v2 + rhs_min. - precedences->Add(v2, v1, MathUtil::CeilOfRatio(rhs_min, magnitude)); - } + LinearExpression2 expr(vars[0], vars[1], coeffs[0], coeffs[1]); + precedences->AddBounds(expr, rhs_min, rhs_max); } else if (vars.size() == 3) { + // TODO(user): This is a weaker duplication of the logic of + // BinaryRelationsMaps, but is is useful for the transitive closure in + // PrecedenceRelations::Build(). Replace this by getting the + // BinaryRelationsMaps affine bounds at level zero. for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { if (i == j) continue; - if (std::abs(coeffs[i]) != std::abs(coeffs[j])) continue; const int other = 3 - i - j; // i + j + other = 0 + 1 + 2. - // Make the terms magnitude * v1 - magnitude * v2 ... - const int64_t magnitude = std::abs(coeffs[i]); - IntegerVariable v1 = vars[i]; - IntegerVariable v2 = vars[j]; - if (coeffs[i] < 0) v1 = NegationOf(v1); - if (coeffs[j] > 0) v2 = NegationOf(v2); - - // magnitude * v1 + other_lb <= magnitude * v2 + rhs_max const int64_t coeff = coeffs[other]; const int64_t other_lb = coeff > 0 ? coeff * integer_trail->LowerBound(vars[other]).value() : coeff * integer_trail->UpperBound(vars[other]).value(); - precedences->Add( - v1, v2, MathUtil::CeilOfRatio(other_lb - rhs_max, magnitude)); - - // magnitude * v1 + other_ub >= magnitude * v2 + rhs_min const int64_t other_ub = coeff > 0 ? coeff * integer_trail->UpperBound(vars[other]).value() : coeff * integer_trail->LowerBound(vars[other]).value(); - precedences->Add( - v2, v1, MathUtil::CeilOfRatio(rhs_min - other_ub, magnitude)); + LinearExpression2 expr(vars[i], vars[j], coeffs[i], coeffs[j]); + precedences->AddBounds(expr, rhs_min - other_ub, rhs_max - other_lb); } } } diff --git a/ortools/sat/cp_model_postsolve.cc b/ortools/sat/cp_model_postsolve.cc index 786e43f660..8231dc2ef1 100644 --- a/ortools/sat/cp_model_postsolve.cc +++ b/ortools/sat/cp_model_postsolve.cc @@ -49,8 +49,12 @@ void PostsolveClause(const ConstraintProto& ct, std::vector* domains) { satisfied = true; } } else { - // We still need to assign free variable. Any value should work. - (*domains)[PositiveRef(ref)] = Domain(0); + // We still need to assign free variable. + // + // It is important to set its value so that the literal in the clause is + // false, so that we support the "filter_sat_postsolve_clauses" option and + // we use a bit less memory for postsolve clauses. + (*domains)[PositiveRef(ref)] = Domain(RefIsPositive(ref) ? 0 : 1); } } if (satisfied) return; diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 76db43e093..d835856bf1 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -2564,7 +2564,7 @@ bool CpModelPresolver::RemoveSingletonInLinear(ConstraintProto* ct) { // also force the postsolve to take search decisions... if (absl::GetFlag(FLAGS_cp_model_debug_postsolve)) { auto* new_ct = context_->NewMappingConstraint(*ct, __FILE__, __LINE__); - const std::string name = new_ct->name(); + const std::string name(new_ct->name()); *new_ct = *ct; new_ct->set_name(absl::StrCat(ct->name(), " copy ", name)); } else { @@ -8221,10 +8221,15 @@ bool CpModelPresolver::PresolvePureSatPart() { // We also disable this if the user asked for tightened domain as this might // fix variable to a potentially infeasible value, and just correct them later // during postsolve of a particular solution. - SatParameters params = context_->params(); - if (params.debug_postsolve_with_full_solver() || - params.fill_tightened_domains_in_response()) { - params.set_presolve_blocked_clause(false); + SatParameters sat_params = context_->params(); + if (sat_params.debug_postsolve_with_full_solver() || + sat_params.fill_tightened_domains_in_response()) { + sat_params.set_presolve_blocked_clause(false); + } + + // This option is only supported by the custom postsolve code. + if (!sat_params.debug_postsolve_with_full_solver()) { + sat_params.set_filter_sat_postsolve_clauses(true); } SatPostsolver sat_postsolver(num_variables); @@ -8280,7 +8285,7 @@ bool CpModelPresolver::PresolvePureSatPart() { // TODO(user): BVA takes time and does not seems to help on the minizinc // benchmarks. So we currently disable it, except if we are on a pure-SAT // problem, where we follow the default (true) or the user specified value. - params.set_presolve_use_bva(false); + sat_params.set_presolve_use_bva(false); } // Disable BVA if we want to keep the symmetries. @@ -8289,7 +8294,7 @@ bool CpModelPresolver::PresolvePureSatPart() { // and also update the generators to take into account the new variables. This // do not seems that easy. if (context_->params().keep_symmetry_in_presolve()) { - params.set_presolve_use_bva(false); + sat_params.set_presolve_use_bva(false); } // Update the time limit of the initial propagation. @@ -8304,7 +8309,7 @@ bool CpModelPresolver::PresolvePureSatPart() { sat_presolver.SetEquivalentLiteralMapping(equiv_map); } sat_presolver.SetTimeLimit(time_limit_); - sat_presolver.SetParameters(params); + sat_presolver.SetParameters(sat_params); // Load in the presolver. // Register the fixed variables with the postsolver. @@ -8354,12 +8359,6 @@ bool CpModelPresolver::PresolvePureSatPart() { for (int i = 0; i < num_variables; ++i) { const int var = new_to_old_index[i]; if (context_->VarToConstraints(var).empty()) { - // Such variable needs to be fixed to some value for the SAT postsolve to - // work. - if (!context_->IsFixed(var)) { - CHECK(context_->IntersectDomainWith( - var, Domain(context_->DomainOf(var).SmallestValue()))); - } context_->MarkVariableAsRemoved(var); } } diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 2928764a79..57a9fc7bcd 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -2282,7 +2282,7 @@ void ParseFromStringOrDie(absl::string_view str, T* proto) { // TODO(user): Support it on android. std::function NewSatParameters( - const std::string& params) { + absl::string_view params) { sat::SatParameters parameters; if (!params.empty()) { ParseFromStringOrDie(params, ¶meters); diff --git a/ortools/sat/cp_model_solver.h b/ortools/sat/cp_model_solver.h index 2afcac1ec6..79d3f7195b 100644 --- a/ortools/sat/cp_model_solver.h +++ b/ortools/sat/cp_model_solver.h @@ -18,6 +18,7 @@ #include #include "absl/flags/declare.h" +#include "absl/strings/string_view.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_parameters.pb.h" @@ -123,8 +124,7 @@ std::function NewBestBoundCallback( \endcode * before calling \c SolveCpModel(). */ -std::function NewSatParameters( - const std::string& params); +std::function NewSatParameters(absl::string_view params); std::function NewSatParameters( const SatParameters& parameters); diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index 23e85a04ff..ab5b398d19 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -33,6 +33,7 @@ #include "absl/numeric/bits.h" #include "absl/types/span.h" // #include "ortools/base/stl_util.h" +#include "ortools/sat/2d_distances_propagator.h" #include "ortools/sat/2d_mandatory_overlap_propagator.h" #include "ortools/sat/2d_orthogonal_packing.h" #include "ortools/sat/2d_try_edge_propagator.h" @@ -187,6 +188,13 @@ void AddNonOverlappingRectangles(const std::vector& x, model->GetOrCreate(); CreateAndRegisterMandatoryOverlapPropagator(helper_2d, model, watcher, 3); + if (model->GetOrCreate() + ->use_linear3_for_no_overlap_2d_precedences()) { + Precedences2DPropagator* propagator = + new Precedences2DPropagator(helper_2d, model); + watcher->SetPropagatorPriority(propagator->RegisterWith(watcher), 4); + model->TakeOwnership(propagator); + } NonOverlappingRectanglesDisjunctivePropagator* constraint = new NonOverlappingRectanglesDisjunctivePropagator(helper_2d, model); diff --git a/ortools/sat/diffn_cuts.cc b/ortools/sat/diffn_cuts.cc index 3c514bd78b..91485e84d9 100644 --- a/ortools/sat/diffn_cuts.cc +++ b/ortools/sat/diffn_cuts.cc @@ -668,11 +668,11 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( if (!helper->SynchronizeAndSetDirection(false, false, false)) { return false; } - generate_cuts("NoOverlap2dXCompletionTime_mirror"); + generate_cuts("NoOverlap2dXCompletionTime"); if (!helper->SynchronizeAndSetDirection(false, false, true)) { return false; } - generate_cuts("NoOverlap2dYCompletionTime_mirror"); + generate_cuts("NoOverlap2dYCompletionTime"); } return true; }; diff --git a/ortools/sat/disjunctive.cc b/ortools/sat/disjunctive.cc index 256673d15c..0faa84e760 100644 --- a/ortools/sat/disjunctive.cc +++ b/ortools/sat/disjunctive.cc @@ -1259,8 +1259,8 @@ bool DisjunctivePrecedences::PropagateSubwindow() { // the offset as much as possible. Note that the alternative of storing it // in PrecedenceData is not necessarily better and harder to update as we // dive/backtrack. - const IntegerValue inner_offset = - precedence_relations_->GetConditionalOffset(end_exp.var, var); + const IntegerValue inner_offset = -precedence_relations_->UpperBound( + LinearExpression2::Difference(end_exp.var, var)); DCHECK_NE(inner_offset, kMinIntegerValue); // We have var >= end_exp.var + inner_offset, so @@ -1312,11 +1312,11 @@ bool DisjunctivePrecedences::PropagateSubwindow() { // Fetch the explanation. // This is okay if a bit slow since we only do that when we push. const AffineExpression& end_exp = helper_->Ends()[ct]; - for (const Literal l : - precedence_relations_->GetConditionalEnforcements(end_exp.var, - var)) { - helper_->MutableLiteralReason()->push_back(l.Negated()); - } + const LinearExpression2 expr = + LinearExpression2::Difference(end_exp.var, var); + precedence_relations_->AddReasonForUpperBoundLowerThan( + expr, precedence_relations_->UpperBound(expr), + helper_->MutableLiteralReason(), helper_->MutableIntegerReason()); } ++stats_.num_propagations; if (!helper_->PushIntegerLiteral( diff --git a/ortools/sat/disjunctive_test.cc b/ortools/sat/disjunctive_test.cc index cf8b2f30cb..85c1db3204 100644 --- a/ortools/sat/disjunctive_test.cc +++ b/ortools/sat/disjunctive_test.cc @@ -249,7 +249,8 @@ TEST(DisjunctiveConstraintTest, Precedences) { CHECK_EQ(e2.coeff, 1); precedences->AddPrecedenceWithOffset(e1.var, e2.var, e1.constant - e2.constant); - relations->Add(e1.var, e2.var, e1.constant - e2.constant); + relations->AddUpperBound(LinearExpression2::Difference(e1.var, e2.var), + e2.constant - e1.constant); }; const int kStart(0); diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index ecfafd04c6..1c8cbf1438 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -531,6 +531,10 @@ class IntegerTrail final : public SatPropagator { IntegerValue LevelZeroLowerBound(AffineExpression exp) const; IntegerValue LevelZeroUpperBound(AffineExpression exp) const; + // Returns globally valid lower/upper bound on the given linear expression. + IntegerValue LevelZeroLowerBound(LinearExpression2 expr) const; + IntegerValue LevelZeroUpperBound(LinearExpression2 expr) const; + // Returns true if the variable is fixed at level 0. bool IsFixedAtLevelZero(IntegerVariable var) const; @@ -1432,6 +1436,30 @@ inline IntegerValue IntegerTrail::LevelZeroUpperBound( return expr.ValueAt(LevelZeroUpperBound(expr.var)); } +inline IntegerValue IntegerTrail::LevelZeroLowerBound( + LinearExpression2 expr) const { + expr.SimpleCanonicalization(); + IntegerValue result = 0; + for (int i = 0; i < 2; ++i) { + if (expr.coeffs[i] != 0) { + result += expr.coeffs[i] * LevelZeroLowerBound(expr.vars[i]); + } + } + return result; +} + +inline IntegerValue IntegerTrail::LevelZeroUpperBound( + LinearExpression2 expr) const { + expr.SimpleCanonicalization(); + IntegerValue result = 0; + for (int i = 0; i < 2; ++i) { + if (expr.coeffs[i] != 0) { + result += expr.coeffs[i] * LevelZeroUpperBound(expr.vars[i]); + } + } + return result; +} + inline bool IntegerTrail::IsFixedAtLevelZero(AffineExpression expr) const { if (expr.var == kNoIntegerVariable) return true; return IsFixedAtLevelZero(expr.var); diff --git a/ortools/sat/integer_base.cc b/ortools/sat/integer_base.cc index 4f3f7e70e7..8af3a695ca 100644 --- a/ortools/sat/integer_base.cc +++ b/ortools/sat/integer_base.cc @@ -56,6 +56,29 @@ void LinearExpression2::SimpleCanonicalization() { } } +IntegerValue LinearExpression2::DivideByGcd() { + const uint64_t gcd = std::gcd(coeffs[0].value(), coeffs[1].value()); + if (gcd > 1) { + coeffs[0] /= gcd; + coeffs[1] /= gcd; + return IntegerValue(gcd); + } + return IntegerValue(1); +} + +bool LinearExpression2::NegateForCanonicalization() { + bool negate = false; + if (coeffs[0] == 0) { + if (coeffs[1] != 0) { + negate = !VariableIsPositive(vars[1]); + } + } else { + negate = !VariableIsPositive(vars[0]); + } + if (negate) Negate(); + return negate; +} + void LinearExpression2::CanonicalizeAndUpdateBounds(IntegerValue& lb, IntegerValue& ub, bool allow_negation) { @@ -63,17 +86,8 @@ void LinearExpression2::CanonicalizeAndUpdateBounds(IntegerValue& lb, if (coeffs[0] == 0 || coeffs[1] == 0) return; // abort. if (allow_negation) { - bool negate = false; - if (coeffs[0] == 0) { - if (coeffs[1] != 0) { - negate = !VariableIsPositive(vars[1]); - } - } else { - negate = !VariableIsPositive(vars[0]); - } - if (negate) { - Negate(); - + const bool negated = NegateForCanonicalization(); + if (negated) { // We need to be able to negate without overflow. CHECK_GE(lb, kMinIntegerValue); CHECK_LE(ub, kMaxIntegerValue); @@ -134,4 +148,21 @@ RelationStatus BestBinaryRelationBounds::GetStatus(LinearExpression2 expr, return RelationStatus::IS_UNKNOWN; } +IntegerValue BestBinaryRelationBounds::GetUpperBound( + LinearExpression2 expr) const { + expr.SimpleCanonicalization(); + const IntegerValue gcd = expr.DivideByGcd(); + const bool negated = expr.NegateForCanonicalization(); + const auto it = best_bounds_.find(expr); + if (it != best_bounds_.end()) { + const auto [known_lb, known_ub] = it->second; + if (negated) { + return CapProdI(gcd, -known_lb); + } else { + return CapProdI(gcd, known_ub); + } + } + return kMaxIntegerValue; +} + } // namespace operations_research::sat diff --git a/ortools/sat/integer_base.h b/ortools/sat/integer_base.h index 3ae647f744..a86d15eb07 100644 --- a/ortools/sat/integer_base.h +++ b/ortools/sat/integer_base.h @@ -346,6 +346,20 @@ H AbslHashValue(H h, const AffineExpression& e) { // A linear expression with at most two variables (coeffs can be zero). // And some utility to canonicalize them. struct LinearExpression2 { + LinearExpression2() = default; + LinearExpression2(IntegerVariable v1, IntegerVariable v2, IntegerValue c1, + IntegerValue c2) { + vars[0] = v1; + vars[1] = v2; + coeffs[0] = c1; + coeffs[1] = c2; + } + + // Build (v1 - v2) + static LinearExpression2 Difference(IntegerVariable v1, IntegerVariable v2) { + return LinearExpression2(v1, v2, 1, -1); + } + // Take the negation of this expression. void Negate() { vars[0] = NegationOf(vars[0]); @@ -355,17 +369,48 @@ struct LinearExpression2 { // This will not change any bounds on the LinearExpression2. // That is we will not potentially Negate() the expression like // CanonicalizeAndUpdateBounds() might do. + // Note that since kNoIntegerVariable=-1 and we sort the variables, if we any + // one zero and one non-zero we will always have the zero first. void SimpleCanonicalization(); - // This fully canonicalize this, and update the given bounds accordingly. + // Fully canonicalizes the expression and updates the given bounds + // accordingly. This is the same as SimpleCanonicalization(), DivideByGcd() + // and the NegateForCanonicalization() with a proper updates of the bounds. void CanonicalizeAndUpdateBounds(IntegerValue& lb, IntegerValue& ub, bool allow_negation = false); + // Divides the expression by the gcd of both coefficients, and returns it. + // Note that we always return something >= 1 even if both coefficients are + // zero. + IntegerValue DivideByGcd(); + + // Makes sure expr and -expr have the same canonical representation by + // negating the expression of it is in the non-canonical form. Returns true if + // the expression was negated. + bool NegateForCanonicalization(); + + absl::Span non_zero_vars() const { + const int first = coeffs[0] == 0 ? 1 : 0; + const int last = coeffs[1] == 0 ? 0 : 1; + return absl::MakeSpan(&vars[first], last - first + 1); + } + + absl::Span non_zero_coeffs() const { + const int first = coeffs[0] == 0 ? 1 : 0; + const int last = coeffs[1] == 0 ? 0 : 1; + return absl::MakeSpan(&coeffs[first], last - first + 1); + } + bool operator==(const LinearExpression2& o) const { return vars[0] == o.vars[0] && vars[1] == o.vars[1] && coeffs[0] == o.coeffs[0] && coeffs[1] == o.coeffs[1]; } + bool operator<(const LinearExpression2& o) const { + return std::tie(vars[0], vars[1], coeffs[0], coeffs[1]) < + std::tie(o.vars[0], o.vars[1], o.coeffs[0], o.coeffs[1]); + } + IntegerValue coeffs[2]; IntegerVariable vars[2]; }; @@ -400,6 +445,11 @@ class BestBinaryRelationBounds { RelationStatus GetStatus(LinearExpression2 expr, IntegerValue lb, IntegerValue ub) const; + // Return a valid upper-bound on the given LinearExpression2. Note that we + // assume kMaxIntegerValue is always valid and returns it if we don't have an + // entry in the hash-map. + IntegerValue GetUpperBound(LinearExpression2 expr) const; + private: // The best bound on the given "canonicalized" expression. absl::flat_hash_map> diff --git a/ortools/sat/integer_base_test.cc b/ortools/sat/integer_base_test.cc index 15f4b987fb..e3b069cd02 100644 --- a/ortools/sat/integer_base_test.cc +++ b/ortools/sat/integer_base_test.cc @@ -37,6 +37,21 @@ TEST(CanonicalizeAffinePrecedenceTest, Basic) { EXPECT_EQ(ub, 5); } +TEST(CanonicalizeAffinePrecedenceTest, OneSingleVariable) { + LinearExpression2 expr; + expr.vars[0] = IntegerVariable(0); + expr.vars[1] = IntegerVariable(0); + expr.coeffs[0] = IntegerValue(2); + expr.coeffs[1] = IntegerValue(2); + + expr.SimpleCanonicalization(); + + EXPECT_EQ(expr.vars[0], kNoIntegerVariable); + EXPECT_EQ(expr.vars[1], IntegerVariable(0)); + EXPECT_EQ(expr.coeffs[0], IntegerValue(0)); + EXPECT_EQ(expr.coeffs[1], IntegerValue(4)); +} + TEST(BestBinaryRelationBoundsTest, Basic) { LinearExpression2 expr; expr.vars[0] = IntegerVariable(0); @@ -63,5 +78,25 @@ TEST(BestBinaryRelationBoundsTest, Basic) { best_bounds.GetStatus(expr, IntegerValue(-5), IntegerValue(3))); } +TEST(BestBinaryRelationBoundsTest, UpperBound) { + LinearExpression2 expr; + expr.vars[0] = IntegerVariable(0); + expr.vars[1] = IntegerVariable(2); + expr.coeffs[0] = IntegerValue(1); + expr.coeffs[1] = IntegerValue(-1); + + BestBinaryRelationBounds best_bounds; + EXPECT_TRUE(best_bounds.Add(expr, IntegerValue(0), IntegerValue(5))); + + EXPECT_EQ(best_bounds.GetUpperBound(expr), IntegerValue(5)); + + expr.coeffs[0] *= 3; + expr.coeffs[1] *= 3; + EXPECT_EQ(best_bounds.GetUpperBound(expr), IntegerValue(15)); + + expr.Negate(); + EXPECT_EQ(best_bounds.GetUpperBound(expr), IntegerValue(0)); +} + } // namespace } // namespace operations_research::sat diff --git a/ortools/sat/intervals.cc b/ortools/sat/intervals.cc index bc5ad0bde2..c50429e71f 100644 --- a/ortools/sat/intervals.cc +++ b/ortools/sat/intervals.cc @@ -155,9 +155,9 @@ IntervalsRepository::GetOrCreateDisjunctivePrecedenceLiteralIfNonTrivial( } // Abort if the relation is already known. - if (relations_maps_->GetPrecedenceStatus(a.end, b.start) == + if (relations_maps_->GetLevelZeroPrecedenceStatus(a.end, b.start) == RelationStatus::IS_TRUE || - relations_maps_->GetPrecedenceStatus(b.end, a.start) == + relations_maps_->GetLevelZeroPrecedenceStatus(b.end, a.start) == RelationStatus::IS_TRUE) { return kNoLiteralIndex; } @@ -217,7 +217,7 @@ bool IntervalsRepository::CreatePrecedenceLiteralIfNonTrivial( // We want l => x <= y and not(l) => x > y <=> y + 1 <= x // Do not create l if the relation is always true or false. - if (relations_maps_->GetPrecedenceStatus(x, y) != + if (relations_maps_->GetLevelZeroPrecedenceStatus(x, y) != RelationStatus::IS_UNKNOWN) { return false; } diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index fde18d9794..fa2b681949 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -1851,7 +1851,8 @@ void LinearProgrammingConstraint::AddMirCuts() { // But there is some degenerate problem where these rows have a really low // weight (or even zero), and having only weight of exactly zero in // std::discrete_distribution will result in a crash. - row_weights[row] = std::max(1e-8, std::abs(simplex_.GetDualValue(row))); + row_weights[row] = + std::max(Fractional(1e-8), std::abs(simplex_.GetDualValue(row))); } // The code here can be really slow, so we put a limit on the number of diff --git a/ortools/sat/linear_propagation.cc b/ortools/sat/linear_propagation.cc index 7632da7836..31fb792911 100644 --- a/ortools/sat/linear_propagation.cc +++ b/ortools/sat/linear_propagation.cc @@ -385,6 +385,7 @@ LinearPropagator::LinearPropagator(Model* model) rev_integer_value_repository_( model->GetOrCreate()), precedences_(model->GetOrCreate()), + binary_relations_(model->GetOrCreate()), random_(model->GetOrCreate()), shared_stats_(model->GetOrCreate()), watcher_id_(watcher_->Register(this)), @@ -534,6 +535,7 @@ bool LinearPropagator::Propagate() { // - Z + Y >= 6 ==> Z >= 1 // - (1) again to push T <= 10 and reach the propagation fixed point. Bitset64::View in_queue = in_queue_.view(); + const bool push_affine_ub = push_affine_ub_for_binary_relations_; while (true) { // We always process the whole queue in FIFO order. // Note that the order really only matter for infeasible constraint so it @@ -571,6 +573,58 @@ bool LinearPropagator::Propagate() { order_.Register(id, NegationOf(var), -new_ub); } } + + // Look at linear3 and update our "linear2 affine upper bound". If we are + // here it means the constraint was in the queue, and its slack changed, + // so it might lead to stronger affine ub. + // + // TODO(user): This can be costly for no reason if we keep updating the + // bound for variable appearing in a single linear3. On another hand it is + // O(1) compared to what this class already do. Profile will tell if it is + // worth it. Maybe we can only share LinearExpression2 that we might look + // up. + // + // TODO(user): This only look at non-enforced linear3. We could look at + // constraint whose enforcement or other variables are fixed at level + // zero, but it is trickier. It could be done if we add a "batch clean up" + // to this class that runs at level zero, and reduce constraints + // accordingly. + const ConstraintInfo& info = infos_[id]; + if (push_affine_ub && info.initial_size == 3 && info.enf_id == -1) { + // A constraint A + B + C <= rhs can lead to up to 3 relations... + const auto vars = GetVariables(info); + const auto coeffs = GetCoeffs(info); + + // We don't "push" relation A + B <= ub if A or B is fixed, because + // the variable bound of the non-fixed A or B should just be as-strong + // as what can be inferred from the binary relation. + if (info.rev_size == 2) { + LinearExpression2 expr; + expr.vars[0] = vars[0]; + expr.vars[1] = vars[1]; + expr.coeffs[0] = coeffs[0]; + expr.coeffs[1] = coeffs[1]; + + // The fixed variable is always at index 2. + // The rev_rhs was updated to: initial_rhs - lb(vars[2]) * coeffs[2]. + const IntegerValue initial_rhs = + info.rev_rhs + coeffs[2] * integer_trail_->LowerBound(vars[2]); + binary_relations_->AddAffineUpperBound( + expr, AffineExpression(vars[2], -coeffs[2], initial_rhs)); + } else if (info.rev_size == 3) { + for (int i = 0; i < 3; ++i) { + LinearExpression2 expr; + const int a = (i + 1) % 3; + const int b = (i + 2) % 3; + expr.vars[0] = vars[a]; + expr.vars[1] = vars[b]; + expr.coeffs[0] = coeffs[a]; + expr.coeffs[1] = coeffs[b]; + binary_relations_->AddAffineUpperBound( + expr, AffineExpression(vars[i], -coeffs[i], info.rev_rhs)); + } + } + } } const int next_id = order_.NextId(); @@ -645,7 +699,7 @@ bool LinearPropagator::AddConstraint( } // Initialize watchers. - // Initialy we want everything to be propagated at least once. + // Initially we want everything to be propagated at least once. in_queue_.resize(in_queue_.size() + 1); if (!enforcement_literals.empty()) { @@ -670,11 +724,13 @@ bool LinearPropagator::AddConstraint( // variables. if (status == EnforcementStatus::IS_ENFORCED) { const auto info = infos_[id]; - if (info.initial_size == 2 && info.all_coeffs_are_one) { + if (info.initial_size == 2) { const auto vars = GetVariables(info); + const auto coeffs = GetCoeffs(info); precedences_->PushConditionalRelation( enforcement_propagator_->GetEnforcementLiterals(enf_id), - vars[0], vars[1], initial_rhs_[id]); + LinearExpression2(vars[0], vars[1], coeffs[0], coeffs[1]), + initial_rhs_[id]); } } }); diff --git a/ortools/sat/linear_propagation.h b/ortools/sat/linear_propagation.h index fd7d375e2d..561ecdf268 100644 --- a/ortools/sat/linear_propagation.h +++ b/ortools/sat/linear_propagation.h @@ -324,6 +324,10 @@ class LinearPropagator : public PropagatorInterface, std::vector* literals_reason, std::vector* trail_indices_reason) final; + void SetPushAffineUbForBinaryRelation() { + push_affine_ub_for_binary_relations_ = true; + } + private: // We try to pack the struct as much as possible. Using a maximum size of // 1 << 29 should be okay since we split long constraint anyway. Technically @@ -389,10 +393,13 @@ class LinearPropagator : public PropagatorInterface, RevIntRepository* rev_int_repository_; RevIntegerValueRepository* rev_integer_value_repository_; PrecedenceRelations* precedences_; + BinaryRelationsMaps* binary_relations_; ModelRandomGenerator* random_; SharedStatistics* shared_stats_ = nullptr; const int watcher_id_; + bool push_affine_ub_for_binary_relations_ = false; + // To know when we backtracked. See SetLevel(). int previous_level_ = 0; diff --git a/ortools/sat/lp_utils.cc b/ortools/sat/lp_utils.cc index a4c84c023d..caf9478faa 100644 --- a/ortools/sat/lp_utils.cc +++ b/ortools/sat/lp_utils.cc @@ -23,6 +23,7 @@ #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/logging.h" #include "ortools/base/strong_vector.h" @@ -1052,7 +1053,7 @@ bool ConvertMPModelProtoToCpModelProto(const SatParameters& params, } case MPGeneralConstraintProto::kAndConstraint: { const auto& and_constraint = general_constraint.and_constraint(); - const std::string& name = general_constraint.name(); + absl::string_view name = general_constraint.name(); ConstraintProto* ct_pos = cp_model->add_constraints(); ct_pos->set_name(name.empty() ? "" : absl::StrCat(name, "_pos")); @@ -1071,7 +1072,7 @@ bool ConvertMPModelProtoToCpModelProto(const SatParameters& params, } case MPGeneralConstraintProto::kOrConstraint: { const auto& or_constraint = general_constraint.or_constraint(); - const std::string& name = general_constraint.name(); + absl::string_view name = general_constraint.name(); ConstraintProto* ct_pos = cp_model->add_constraints(); ct_pos->set_name(name.empty() ? "" : absl::StrCat(name, "_pos")); diff --git a/ortools/sat/precedences.cc b/ortools/sat/precedences.cc index 0b6c0e0ed6..562795dd3c 100644 --- a/ortools/sat/precedences.cc +++ b/ortools/sat/precedences.cc @@ -53,38 +53,60 @@ namespace operations_research { namespace sat { -bool PrecedenceRelations::Add(IntegerVariable tail, IntegerVariable head, - IntegerValue offset) { - // Ignore trivial relation: tail + offset <= head. - if (integer_trail_->LevelZeroUpperBound(tail) + offset <= - integer_trail_->LevelZeroLowerBound(head)) { +bool PrecedenceRelations::AddBounds(LinearExpression2 expr, IntegerValue lb, + IntegerValue ub) { + expr.CanonicalizeAndUpdateBounds(lb, ub); + + if (expr.coeffs[0] == 0 || expr.coeffs[1] == 0) { + // This class handles only binary relationships, let something else handle + // the case where there is actually a single variable. return false; } - // TODO(user): Return infeasible if tail == head and offset > 0. - // TODO(user): if tail = Negation(head) also update Domain. - if (tail == head) return false; - // Add to root_relations_. // // TODO(user): AddInternal() only returns true if this is the first relation // between head and tail. But we can still avoid an extra lookup. - if (offset <= GetOffset(tail, head)) return false; - AddInternal(tail, head, offset); + const bool add_ub = ub < LevelZeroUpperBound(expr); + LinearExpression2 expr_for_lb = expr; + expr_for_lb.Negate(); + const bool add_lb = lb > -LevelZeroUpperBound(expr_for_lb); + if (!add_ub && !add_lb) { + return false; + } + + if (add_ub) { + AddInternal(expr, ub); + } + if (add_lb) { + AddInternal(expr_for_lb, -lb); + } // If we are not built, make sure there is enough room in the graph. // TODO(user): Alternatively, force caller to do a Resize(). const int max_node = - std::max(PositiveVariable(tail), PositiveVariable(head)).value() + 1; + std::max(PositiveVariable(expr.vars[0]), PositiveVariable(expr.vars[1])) + .value() + + 1; if (!is_built_ && max_node >= graph_.num_nodes()) { graph_.AddNode(max_node); } return true; } +bool PrecedenceRelations::AddUpperBound(LinearExpression2 expr, + IntegerValue ub) { + return AddBounds(expr, kMinIntegerValue, ub); +} + void PrecedenceRelations::PushConditionalRelation( - absl::Span enforcements, IntegerVariable a, - IntegerVariable b, IntegerValue rhs) { + absl::Span enforcements, LinearExpression2 expr, + IntegerValue rhs) { + expr.SimpleCanonicalization(); + if (expr.coeffs[0] == 0 || expr.coeffs[1] == 0) { + return; + } + // This must be currently true. if (DEBUG_MODE) { for (const Literal l : enforcements) { @@ -93,14 +115,16 @@ void PrecedenceRelations::PushConditionalRelation( } if (enforcements.empty() || trail_->CurrentDecisionLevel() == 0) { - Add(a, NegationOf(b), -rhs); + AddUpperBound(expr, rhs); return; } + const IntegerValue gcd = expr.DivideByGcd(); + rhs = FloorRatio(rhs, gcd); + // Ignore if no better than best_relations, otherwise increase it. - const auto key = GetKey(a, b); { - const auto [it, inserted] = best_relations_.insert({key, rhs}); + const auto [it, inserted] = best_relations_.insert({expr, rhs}); if (!inserted) { if (rhs >= it->second) return; // Ignore. it->second = rhs; @@ -108,17 +132,20 @@ void PrecedenceRelations::PushConditionalRelation( } const int new_index = conditional_stack_.size(); - const auto [it, inserted] = conditional_relations_.insert({key, new_index}); + const auto [it, inserted] = conditional_relations_.insert({expr, new_index}); if (inserted) { CreateLevelEntryIfNeeded(); - conditional_stack_.emplace_back(/*prev_entry=*/-1, rhs, key, enforcements); + conditional_stack_.emplace_back(/*prev_entry=*/-1, rhs, expr, enforcements); - const int new_size = std::max(a.value(), b.value()) + 1; - if (new_size > conditional_after_.size()) { - conditional_after_.resize(new_size); + if (expr.coeffs[0] == 1 && expr.coeffs[1] == 1) { + const int new_size = + std::max(expr.vars[0].value(), expr.vars[1].value()) + 1; + if (new_size > conditional_after_.size()) { + conditional_after_.resize(new_size); + } + conditional_after_[expr.vars[0]].push_back(NegationOf(expr.vars[1])); + conditional_after_[expr.vars[1]].push_back(NegationOf(expr.vars[0])); } - conditional_after_[a].push_back(NegationOf(b)); - conditional_after_[b].push_back(NegationOf(a)); } else { // We should only decrease because we ignored entry worse than the one in // best_relations_. @@ -128,7 +155,7 @@ void PrecedenceRelations::PushConditionalRelation( // Update. it->second = new_index; CreateLevelEntryIfNeeded(); - conditional_stack_.emplace_back(prev_entry, rhs, key, enforcements); + conditional_stack_.emplace_back(prev_entry, rhs, expr, enforcements); } } @@ -155,12 +182,14 @@ void PrecedenceRelations::SetLevel(int level) { UpdateBestRelation(back.key, kMaxIntegerValue); conditional_relations_.erase(back.key); - DCHECK_EQ(conditional_after_[back.key.first].back(), - NegationOf(back.key.second)); - DCHECK_EQ(conditional_after_[back.key.second].back(), - NegationOf(back.key.first)); - conditional_after_[back.key.first].pop_back(); - conditional_after_[back.key.second].pop_back(); + if (back.key.coeffs[0] == 1 && back.key.coeffs[1] == 1) { + DCHECK_EQ(conditional_after_[back.key.vars[0]].back(), + NegationOf(back.key.vars[1])); + DCHECK_EQ(conditional_after_[back.key.vars[1]].back(), + NegationOf(back.key.vars[0])); + conditional_after_[back.key.vars[0]].pop_back(); + conditional_after_[back.key.vars[1]].pop_back(); + } } conditional_stack_.pop_back(); } @@ -168,19 +197,26 @@ void PrecedenceRelations::SetLevel(int level) { } } -IntegerValue PrecedenceRelations::GetOffset(IntegerVariable a, - IntegerVariable b) const { - const auto it = root_relations_.find(GetKey(a, NegationOf(b))); +IntegerValue PrecedenceRelations::LevelZeroUpperBound( + LinearExpression2 expr) const { + expr.SimpleCanonicalization(); + const IntegerValue gcd = expr.DivideByGcd(); + const auto it = root_relations_.find(expr); if (it != root_relations_.end()) { - return -it->second; + return CapProdI(it->second, gcd); } - return kMinIntegerValue; + return kMaxIntegerValue; } -absl::Span PrecedenceRelations::GetConditionalEnforcements( - IntegerVariable a, IntegerVariable b) const { - const auto it = conditional_relations_.find(GetKey(a, NegationOf(b))); - if (it == conditional_relations_.end()) return {}; +void PrecedenceRelations::AddReasonForUpperBoundLowerThan( + LinearExpression2 expr, IntegerValue ub, + std::vector* literal_reason, + std::vector* /*unused*/) const { + expr.SimpleCanonicalization(); + if (ub >= LevelZeroUpperBound(expr)) return; + const IntegerValue gcd = expr.DivideByGcd(); + const auto it = conditional_relations_.find(expr); + DCHECK(it != conditional_relations_.end()); const ConditionalEntry& entry = conditional_stack_[it->second]; if (DEBUG_MODE) { @@ -188,23 +224,23 @@ absl::Span PrecedenceRelations::GetConditionalEnforcements( CHECK(trail_->Assignment().LiteralIsTrue(l)); } } - const IntegerValue root_level_offset = GetOffset(a, b); - const IntegerValue conditional_offset = -entry.rhs; - if (conditional_offset <= root_level_offset) return {}; - - DCHECK_EQ(entry.rhs, -GetConditionalOffset(a, b)); - return entry.enforcements; + DCHECK_EQ(CapProdI(gcd, entry.rhs), UpperBound(expr)); + DCHECK_LE(CapProdI(gcd, entry.rhs), ub); + for (const Literal l : entry.enforcements) { + literal_reason->push_back(l.Negated()); + } } -IntegerValue PrecedenceRelations::GetConditionalOffset( - IntegerVariable a, IntegerVariable b) const { - const auto it = best_relations_.find(GetKey(a, NegationOf(b))); +IntegerValue PrecedenceRelations::UpperBound(LinearExpression2 expr) const { + expr.SimpleCanonicalization(); + const IntegerValue gcd = expr.DivideByGcd(); + const auto it = best_relations_.find(expr); if (it != best_relations_.end()) { - return -it->second; + return CapProdI(gcd, it->second); } - DCHECK(!root_relations_.contains(GetKey(a, NegationOf(b)))); - DCHECK(!conditional_relations_.contains(GetKey(a, NegationOf(b)))); - return kMinIntegerValue; + DCHECK(!root_relations_.contains(expr)); + DCHECK(!conditional_relations_.contains(expr)); + return kMaxIntegerValue; } void PrecedenceRelations::Build() { @@ -219,9 +255,8 @@ void PrecedenceRelations::Build() { // And use this to compute the "closure". CHECK(arc_offsets_.empty()); graph_.ReserveArcs(2 * root_relations_.size()); - std::vector< - std::pair, IntegerValue>> - root_relations_sorted(root_relations_.begin(), root_relations_.end()); + std::vector> root_relations_sorted( + root_relations_.begin(), root_relations_.end()); std::sort(root_relations_sorted.begin(), root_relations_sorted.end()); for (const auto [var_pair, negated_offset] : root_relations_sorted) { // TODO(user): Support negative offset? @@ -232,21 +267,26 @@ void PrecedenceRelations::Build() { const IntegerValue offset = -negated_offset; if (offset < 0) continue; + if (var_pair.coeffs[0] != 1 || var_pair.coeffs[1] != 1) { + // TODO(user): Support non-1 coefficients. + continue; + } + // We have two arcs. { - const IntegerVariable tail = var_pair.first; - const IntegerVariable head = NegationOf(var_pair.second); + const IntegerVariable tail = var_pair.vars[0]; + const IntegerVariable head = NegationOf(var_pair.vars[1]); graph_.AddArc(tail.value(), head.value()); arc_offsets_.push_back(offset); - CHECK_LT(var_pair.second, before.size()); + CHECK_LT(var_pair.vars[1], before.size()); before[head].push_back(tail); } { - const IntegerVariable tail = var_pair.second; - const IntegerVariable head = NegationOf(var_pair.first); + const IntegerVariable tail = var_pair.vars[1]; + const IntegerVariable head = NegationOf(var_pair.vars[0]); graph_.AddArc(tail.value(), head.value()); arc_offsets_.push_back(offset); - CHECK_LT(var_pair.second, before.size()); + CHECK_LT(var_pair.vars[1], before.size()); before[head].push_back(tail); } } @@ -294,16 +334,19 @@ void PrecedenceRelations::Build() { const IntegerValue arc_offset = arc_offsets_[arc]; if (++work > kWorkLimit) break; - if (AddInternal(tail_var, head_var, arc_offset)) { + if (AddInternal(LinearExpression2::Difference(tail_var, head_var), + -arc_offset)) { before[head_var].push_back(tail_var); } for (const IntegerVariable before_var : before[tail_var]) { if (++work > kWorkLimit) break; + LinearExpression2 expr_for_key(before_var, tail_var, 1, -1); + expr_for_key.SimpleCanonicalization(); const IntegerValue offset = - -root_relations_.at(GetKey(before_var, NegationOf(tail_var))) + - arc_offset; - if (AddInternal(before_var, head_var, offset)) { + -root_relations_.at(expr_for_key) + arc_offset; + if (AddInternal(LinearExpression2::Difference(before_var, head_var), + -offset)) { before[head_var].push_back(before_var); } } @@ -582,8 +625,9 @@ void PrecedencesPropagator::PushConditionalRelations(const ArcInfo& arc) { // add this to the reason though. if (arc.offset_var != kNoIntegerVariable) return; const IntegerValue offset = ArcOffset(arc); - relations_->PushConditionalRelation(arc.presence_literals, arc.tail_var, - NegationOf(arc.head_var), -offset); + relations_->PushConditionalRelation( + arc.presence_literals, + LinearExpression2::Difference(arc.tail_var, arc.head_var), -offset); } void PrecedencesPropagator::Untrail(const Trail& trail, int trail_index) { @@ -1490,6 +1534,7 @@ int GreaterThanAtLeastOneOfDetector::AddGreaterThanAtLeastOneOfConstraints( BinaryRelationsMaps::BinaryRelationsMaps(Model* model) : integer_trail_(model->GetOrCreate()), integer_encoder_(model->GetOrCreate()), + watcher_(model->GetOrCreate()), shared_stats_(model->GetOrCreate()) { int index = 0; model->GetOrCreate()->callbacks.push_back( @@ -1524,9 +1569,24 @@ BinaryRelationsMaps::~BinaryRelationsMaps() { if (!VLOG_IS_ON(1)) return; std::vector> stats; stats.push_back({"BinaryRelationsMaps/num_relations", num_updates_}); + stats.push_back( + {"BinaryRelationsMaps/num_affine_updates", num_affine_updates_}); shared_stats_->AddStats(stats); } +IntegerValue BinaryRelationsMaps::GetImpliedUpperBound( + const LinearExpression2& expr) const { + DCHECK_GE(expr.coeffs[0], 0); + DCHECK_GE(expr.coeffs[1], 0); + IntegerValue implied_ub = 0; + for (const int i : {0, 1}) { + if (expr.coeffs[i] > 0) { + implied_ub += expr.coeffs[i] * integer_trail_->UpperBound(expr.vars[i]); + } + } + return implied_ub; +} + std::pair BinaryRelationsMaps::GetImpliedLevelZeroBounds( const LinearExpression2& expr) const { @@ -1561,16 +1621,16 @@ void BinaryRelationsMaps::AddRelationBounds(LinearExpression2 expr, if (lb > ub) return; // unsat ?? if (lb == implied_lb && ub == implied_ub) return; // trivially true. - if (best_upper_bounds_.Add(expr, lb, ub)) { + if (best_root_level_bounds_.Add(expr, lb, ub)) { // TODO(user): Also push them to a global shared repository after // remapping IntegerVariable to proto indices. ++num_updates_; } } -RelationStatus BinaryRelationsMaps::GetStatus(LinearExpression2 expr, - IntegerValue lb, - IntegerValue ub) const { +RelationStatus BinaryRelationsMaps::GetLevelZeroStatus(LinearExpression2 expr, + IntegerValue lb, + IntegerValue ub) const { expr.CanonicalizeAndUpdateBounds(lb, ub); const auto [implied_lb, implied_ub] = GetImpliedLevelZeroBounds(expr); lb = std::max(lb, implied_lb); @@ -1580,11 +1640,11 @@ RelationStatus BinaryRelationsMaps::GetStatus(LinearExpression2 expr, if (lb > ub) return RelationStatus::IS_FALSE; if (lb == implied_lb && ub == implied_ub) return RelationStatus::IS_TRUE; - // Relax as best_upper_bounds_.GetStatus() might have older bounds. + // Relax as best_root_level_bounds_.GetStatus() might have older bounds. if (lb == implied_lb) lb = kMinIntegerValue; if (ub == implied_ub) ub = kMaxIntegerValue; - return best_upper_bounds_.GetStatus(expr, lb, ub); + return best_root_level_bounds_.GetStatus(expr, lb, ub); } std::pair BinaryRelationsMaps::FromDifference( @@ -1600,17 +1660,17 @@ std::pair BinaryRelationsMaps::FromDifference( return {std::move(expr), ub}; } -RelationStatus BinaryRelationsMaps::GetPrecedenceStatus( +RelationStatus BinaryRelationsMaps::GetLevelZeroPrecedenceStatus( AffineExpression a, AffineExpression b) const { const auto [expr, ub] = FromDifference(a, b); - return GetStatus(expr, kMinIntegerValue, ub); + return GetLevelZeroStatus(expr, kMinIntegerValue, ub); } void BinaryRelationsMaps::AddReifiedPrecedenceIfNonTrivial(Literal l, AffineExpression a, AffineExpression b) { const auto [expr, ub] = FromDifference(a, b); - const RelationStatus status = GetStatus(expr, kMinIntegerValue, ub); + const RelationStatus status = GetLevelZeroStatus(expr, kMinIntegerValue, ub); if (status != RelationStatus::IS_UNKNOWN) return; relation_to_lit_.insert({{expr, ub}, l}); @@ -1622,7 +1682,7 @@ void BinaryRelationsMaps::AddReifiedPrecedenceIfNonTrivial(Literal l, LiteralIndex BinaryRelationsMaps::GetReifiedPrecedence(AffineExpression a, AffineExpression b) { const auto [expr, ub] = FromDifference(a, b); - const RelationStatus status = GetStatus(expr, kMinIntegerValue, ub); + const RelationStatus status = GetLevelZeroStatus(expr, kMinIntegerValue, ub); if (status == RelationStatus::IS_TRUE) { return integer_encoder_->GetTrueLiteral().Index(); } @@ -1635,5 +1695,116 @@ LiteralIndex BinaryRelationsMaps::GetReifiedPrecedence(AffineExpression a, return it->second; } +bool BinaryRelationsMaps::AddAffineUpperBound(LinearExpression2 expr, + AffineExpression affine_ub) { + const IntegerValue new_ub = integer_trail_->UpperBound(affine_ub); + expr.SimpleCanonicalization(); + + // Not better than trivial upper bound. + if (GetImpliedUpperBound(expr) <= new_ub) return false; + + // Not better than the root level upper bound. + if (best_root_level_bounds_.GetUpperBound(expr) <= new_ub) return false; + + const IntegerValue gcd = expr.DivideByGcd(); + + const auto it = best_affine_ub_.find(expr); + if (it != best_affine_ub_.end()) { + const auto [old_affine_ub, old_gcd] = it->second; + // We have an affine bound for this expr in the map. Can be exactly the + // same, a better one or a worse one. + if (old_affine_ub == affine_ub && old_gcd == gcd) { + // The affine bound is already in the map. + NotifyWatchingPropagators(); // The affine bound was updated. + return false; + } + const IntegerValue old_ub = + FloorRatio(integer_trail_->UpperBound(old_affine_ub), old_gcd); + if (old_ub <= new_ub) return false; // old bound is better. + } + + // We have gcd * canonical_expr <= affine_ub, so we do need to store a + // "divisor". + ++num_affine_updates_; + best_affine_ub_[expr] = {affine_ub, gcd}; + NotifyWatchingPropagators(); + return true; +} + +void BinaryRelationsMaps::NotifyWatchingPropagators() const { + for (const int id : propagator_ids_) { + watcher_->CallOnNextPropagate(id); + } +} + +IntegerValue BinaryRelationsMaps::UpperBound(LinearExpression2 expr) const { + expr.SimpleCanonicalization(); + + const IntegerValue trivial_ub = GetImpliedUpperBound(expr); + const IntegerValue root_level_ub = + best_root_level_bounds_.GetUpperBound(expr); + const IntegerValue best_ub = std::min(root_level_ub, trivial_ub); + + const IntegerValue gcd = expr.DivideByGcd(); + const auto it = best_affine_ub_.find(expr); + if (it == best_affine_ub_.end()) { + return best_ub; + } else { + const auto [affine, divisor] = it->second; + const IntegerValue canonical_ub = + FloorRatio(integer_trail_->UpperBound(affine), divisor); + return std::min(best_ub, CapProdI(gcd, canonical_ub)); + } +} + +// TODO(user): If the trivial bound is better, its explanation is different... +void BinaryRelationsMaps::AddReasonForUpperBoundLowerThan( + LinearExpression2 expr, IntegerValue ub, + std::vector* /*literal_reason*/, + std::vector* integer_reason) const { + expr.SimpleCanonicalization(); + + if (expr.coeffs[0] == 0 && expr.coeffs[1] == 0) return; // trivially zero + + // Starts by simple bounds. + if (best_root_level_bounds_.GetUpperBound(expr) <= ub) return; + + // Add explanation if it is a trivial bound. + const IntegerValue implied_ub = GetImpliedUpperBound(expr); + if (implied_ub <= ub) { + const IntegerValue slack = ub - implied_ub; + expr.Negate(); // AppendRelaxedLinearReason() explains a lower bound. + absl::Span vars = expr.non_zero_vars(); + absl::Span coeffs = expr.non_zero_coeffs(); + integer_trail_->AppendRelaxedLinearReason(slack, coeffs, vars, + integer_reason); + return; + } + + // None of the bound above are enough, try the affine one. Note that gcd * + // expr <= ub, is the same as asking why expr <= FloorRatio(ub, gcd). + const IntegerValue gcd = expr.DivideByGcd(); + const auto it = best_affine_ub_.find(expr); + if (it == best_affine_ub_.end()) return; + + // We want the reason for "expr <= ub", that is the reason for + // - "gcd * canonical_expr <= ub" + // - "canonical_expr <= FloorRatio(ub, gcd); + // + // knowing that canonical_expr <= affine_ub / divisor. + const auto [affine, divisor] = it->second; + integer_reason->push_back( + affine.LowerOrEqual(CapProdI(FloorRatio(ub, gcd) + 1, divisor) - 1)); +} + +std::vector +BinaryRelationsMaps::GetAllExpressionsWithAffineBounds() const { + std::vector result; + for (const auto [expr, info] : best_affine_ub_) { + result.push_back(expr); + } + return result; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/precedences.h b/ortools/sat/precedences.h index 8056303e86..4ded3dbc4d 100644 --- a/ortools/sat/precedences.h +++ b/ortools/sat/precedences.h @@ -47,7 +47,7 @@ struct FullIntegerPrecedence { std::vector offsets; }; -// Stores all the precedences relation of the form "tail_X + offset <= head_X" +// Stores all the precedences relation of the form "a*x + b*y <= ub" // that we could extract from the linear constraint of the model. These are // stored in a directed graph. // @@ -68,11 +68,15 @@ class PrecedenceRelations : public ReversibleInterface { graph_.AddNode(num_variables - 1); } - // Add a relation tail + offset <= head. - // Returns true if it was added and is considered "new". - bool Add(IntegerVariable tail, IntegerVariable head, IntegerValue offset); + // Add a relation lb <= expr <= ub. If expr is not a proper linear2 expression + // (e.g. 0*x + y, y + y, y - y) it will be ignored. Returns true if it was + // added and is considered "new". + bool AddBounds(LinearExpression2 expr, IntegerValue lb, IntegerValue ub); - // Adds add relation (enf => a + b <= rhs) that is assumed to be true at + // Same as above, but only for the upper bound. + bool AddUpperBound(LinearExpression2 expr, IntegerValue ub); + + // Adds add relation (enf => expr <= rhs) that is assumed to be true at // the current level. // // It will be automatically reverted via the SetLevel() functions that is @@ -80,9 +84,11 @@ class PrecedenceRelations : public ReversibleInterface { // // This is assumed to be called when a relation becomes true (enforcement are // assigned) and when it becomes false in reverse order (CHECKed). + // + // If expr is not a proper linear2 expression (e.g. 0*x + y, y + y, y - y) it + // will be ignored. void PushConditionalRelation(absl::Span enforcements, - IntegerVariable a, IntegerVariable b, - IntegerValue rhs); + LinearExpression2 expr, IntegerValue rhs); // Called each time we change decision level. void SetLevel(int level) final; @@ -92,6 +98,9 @@ class PrecedenceRelations : public ReversibleInterface { // This currently only works if the precedence relation form a DAG. // If not we will just abort. TODO(user): generalize. // + // For more efficiency, this method ignores all linear2 expressions with any + // coefficient different from 1. + // // TODO(user): Put some work limit in place, as this can be slow. Complexity // is in O(vars.size()) * num_arcs. // @@ -106,7 +115,10 @@ class PrecedenceRelations : public ReversibleInterface { // Returns a set of precedences (var, index) such that var is after // vars[index]. All entries for the same variable will be contiguous and // sorted by index. We only list variable with at least two entries. The - // offset can be retrieved via GetConditionalOffset(vars[index], var). + // offset can be retrieved via UpperBound(vars[index], var). + // + // For more efficiency, this method ignores all linear2 expressions with any + // coefficient different from 1. struct PrecedenceData { IntegerVariable var; int index; @@ -120,23 +132,26 @@ class PrecedenceRelations : public ReversibleInterface { // // Warning: If there are too many, this will NOT contain all relations. // - // Returns kMinIntegerValue if there are none. - // Otherwise a + offset <= b. - IntegerValue GetOffset(IntegerVariable a, IntegerVariable b) const; + // Returns kMaxIntegerValue if there are none, otherwise return an upper bound + // such that expr <= ub. + IntegerValue LevelZeroUpperBound(LinearExpression2 expr) const; - // Returns the minimum distance between a and b, and the reason for it (all - // true). Note that we always check GetOffset() so if it is better, the - // returned literal reason will be empty. + // Returns the maximum value for expr, and the reason for it (all + // true). Note that we always check LevelZeroUpperBound() so if it is better, + // the returned literal reason will be empty. // // We separate the two because usually the reason is only needed when we push, // which happen less often, so we don't mind doing two hash lookups, and we - // really want to optimize the GetConditionalOffset() instead. + // really want to optimize the UpperBound() instead. // // Important: This doesn't contains the transitive closure. // Important: The span is only valid in a narrow scope. - IntegerValue GetConditionalOffset(IntegerVariable a, IntegerVariable b) const; - absl::Span GetConditionalEnforcements(IntegerVariable a, - IntegerVariable b) const; + IntegerValue UpperBound(LinearExpression2 expr) const; + + void AddReasonForUpperBoundLowerThan( + LinearExpression2 expr, IntegerValue ub, + std::vector* literal_reason, + std::vector* integer_reason) const; // The current code requires the internal data to be processed once all // relations are loaded. @@ -147,47 +162,45 @@ class PrecedenceRelations : public ReversibleInterface { private: void CreateLevelEntryIfNeeded(); - std::pair GetKey(IntegerVariable a, - IntegerVariable b) const { - return a <= b ? std::make_pair(a, b) : std::make_pair(b, a); - } - - // tail + offset <= head. - // Which is the same as tail - head <= -offset. - bool AddInternal(IntegerVariable tail, IntegerVariable head, - IntegerValue offset) { - const auto key = GetKey(tail, NegationOf(head)); - const auto [it, inserted] = root_relations_.insert({key, -offset}); - UpdateBestRelationIfBetter(key, -offset); + // expr <= ub. + bool AddInternal(LinearExpression2 expr, IntegerValue ub) { + expr.SimpleCanonicalization(); + if (expr.coeffs[0] == 0 || expr.coeffs[1] == 0) { + return false; + } + const auto [it, inserted] = root_relations_.insert({expr, ub}); + UpdateBestRelationIfBetter(expr, ub); if (inserted) { - const int new_size = std::max(tail.value(), NegationOf(head).value()) + 1; + if (expr.coeffs[0] != 1 || expr.coeffs[1] != 1) { + return true; + } + const int new_size = + std::max(expr.vars[0].value(), expr.vars[1].value()) + 1; if (new_size > after_.size()) after_.resize(new_size); - after_[tail].push_back(head); - after_[NegationOf(head)].push_back(NegationOf(tail)); + after_[expr.vars[0]].push_back(NegationOf(expr.vars[1])); + after_[expr.vars[1]].push_back(NegationOf(expr.vars[0])); return true; } - it->second = std::min(it->second, -offset); + it->second = std::min(it->second, ub); return false; } - void UpdateBestRelationIfBetter( - std::pair key, IntegerValue rhs) { - const auto [it, inserted] = best_relations_.insert({key, rhs}); + void UpdateBestRelationIfBetter(LinearExpression2 expr, IntegerValue rhs) { + const auto [it, inserted] = best_relations_.insert({expr, rhs}); if (!inserted) { it->second = std::min(it->second, rhs); } } - void UpdateBestRelation(std::pair key, - IntegerValue rhs) { - const auto it = root_relations_.find(key); + void UpdateBestRelation(LinearExpression2 expr, IntegerValue rhs) { + const auto it = root_relations_.find(expr); if (it != root_relations_.end()) { rhs = std::min(rhs, it->second); } if (rhs == kMaxIntegerValue) { - best_relations_.erase(key); + best_relations_.erase(expr); } else { - best_relations_[key] = rhs; + best_relations_[expr] = rhs; } } @@ -207,34 +220,31 @@ class PrecedenceRelations : public ReversibleInterface { // TODO(user): this kind of reversible hash_map is already implemented in // other part of the code. Consolidate. struct ConditionalEntry { - ConditionalEntry(int p, IntegerValue r, - std::pair k, + ConditionalEntry(int p, IntegerValue r, LinearExpression2 k, absl::Span e) : prev_entry(p), rhs(r), key(k), enforcements(e.begin(), e.end()) {} int prev_entry; IntegerValue rhs; - std::pair key; + LinearExpression2 key; absl::InlinedVector enforcements; }; std::vector conditional_stack_; std::vector> level_to_stack_size_; - // This is always stored in the form (a + b <= rhs). + // This is always stored in the form (expr <= rhs). // The conditional relations contains indices in the conditional_stack_. - absl::flat_hash_map, IntegerValue> - root_relations_; - absl::flat_hash_map, int> - conditional_relations_; + absl::flat_hash_map root_relations_; + absl::flat_hash_map conditional_relations_; // Contains std::min() of the offset from root_relations_ and // conditional_relations_. - absl::flat_hash_map, IntegerValue> - best_relations_; + absl::flat_hash_map best_relations_; - // Store for each variable x, the variables y that appears in GetOffset(x, y) - // or GetConditionalOffset(x, y). That is the variable that are after x with - // an offset. Note that conditional_after_ is updated on dive/backtrack. + // Store for each variable x, the variables y that appears alongside it in + // LevelZeroUpperBound(expr) or UpperBound(expr). That is the variable + // that are after x with an offset. Note that conditional_after_ is updated on + // dive/backtrack. util_intops::StrongVector> after_; util_intops::StrongVector> @@ -588,12 +598,12 @@ class BinaryRelationsMaps { // relation. void AddRelationBounds(LinearExpression2 expr, IntegerValue lb, IntegerValue ub); - RelationStatus GetStatus(LinearExpression2 expr, IntegerValue lb, - IntegerValue ub) const; + RelationStatus GetLevelZeroStatus(LinearExpression2 expr, IntegerValue lb, + IntegerValue ub) const; // Return the status of a <= b; - RelationStatus GetPrecedenceStatus(AffineExpression a, - AffineExpression b) const; + RelationStatus GetLevelZeroPrecedenceStatus(AffineExpression a, + AffineExpression b) const; // Register the fact that l <=> ( a <= b ). // These are considered equivalence relation. @@ -605,20 +615,48 @@ class BinaryRelationsMaps { // true/false literal if the status is known at level zero. LiteralIndex GetReifiedPrecedence(AffineExpression a, AffineExpression b); + // If the given upper bound evaluate better than the current one we have, this + // will replace it and returns true, otherwise it returns false. + // + // Note that we never store trivial upper bound (using the current variable + // domain). + bool AddAffineUpperBound(LinearExpression2 expr, AffineExpression affine_ub); + + // Returns the best known upper-bound of the given LinearExpression2 at the + // current decision level. If its explanation is needed, it can be queried + // with the second function. + IntegerValue UpperBound(LinearExpression2 expr) const; + void AddReasonForUpperBoundLowerThan( + LinearExpression2 expr, IntegerValue ub, + std::vector* literal_reason, + std::vector* integer_reason) const; + + // Warning, the order will not be deterministic. + std::vector GetAllExpressionsWithAffineBounds() const; + + int NumExpressionsWithAffineBounds() const { return best_affine_ub_.size(); } + + void WatchAllLinearExpressions2(int id) { propagator_ids_.insert(id); } + private: + void NotifyWatchingPropagators() const; + // Return the pair (a - b) <= rhs. std::pair FromDifference( const AffineExpression& a, const AffineExpression& b) const; + IntegerValue GetImpliedUpperBound(const LinearExpression2& expr) const; std::pair GetImpliedLevelZeroBounds( const LinearExpression2& expr) const; IntegerTrail* integer_trail_; IntegerEncoder* integer_encoder_; + GenericLiteralWatcher* watcher_; SharedStatistics* shared_stats_; - BestBinaryRelationBounds best_upper_bounds_; + BestBinaryRelationBounds best_root_level_bounds_; int64_t num_updates_ = 0; + int64_t num_affine_updates_ = 0; // This stores relations l <=> (linear2 <= rhs). absl::flat_hash_map, Literal> @@ -631,6 +669,18 @@ class BinaryRelationsMaps { absl::flat_hash_set variable_appearing_in_reified_relations_; std::vector> all_reified_relations_; + + // This stores linear2 <= AffineExpression / divisor. + // + // Note(user): This is a "cheap way" to not have to deal with backtracking, If + // we have many possible AffineExpression that bounds a LinearExpression2, we + // keep the best one during "search dive" but on backtrack we might have a + // sub-optimal relation. + absl::flat_hash_map> + best_affine_ub_; + + absl::btree_set propagator_ids_; }; // Detects if at least one of a subset of linear of size 2 or 1, touching the @@ -731,7 +781,9 @@ inline std::function LowerOrEqualWithOffset(IntegerVariable a, IntegerVariable b, int64_t offset) { return [=](Model* model) { - model->GetOrCreate()->Add(a, b, IntegerValue(offset)); + LinearExpression2 expr(a, b, 1, -1); + model->GetOrCreate()->AddUpperBound( + expr, IntegerValue(-offset)); model->GetOrCreate()->AddPrecedenceWithOffset( a, b, IntegerValue(offset)); }; @@ -745,8 +797,9 @@ inline std::function AffineCoeffOneLowerOrEqualWithOffset( CHECK_NE(b.var, kNoIntegerVariable); CHECK_EQ(b.coeff, 1); return [=](Model* model) { - model->GetOrCreate()->Add( - a.var, b.var, a.constant - b.constant + offset); + LinearExpression2 expr(a.var, b.var, 1, -1); + model->GetOrCreate()->AddUpperBound( + expr, -a.constant + b.constant - offset); model->GetOrCreate()->AddPrecedenceWithOffset( a.var, b.var, a.constant - b.constant + offset); }; @@ -758,8 +811,9 @@ inline void AddConditionalSum2LowerOrEqual( IntegerVariable b, int64_t ub, Model* model) { // TODO(user): Refactor to be sure we do not miss any level zero relations. if (enforcement_literals.empty()) { - model->GetOrCreate()->Add(a, NegationOf(b), - IntegerValue(-ub)); + LinearExpression2 expr(a, b, 1, 1); + model->GetOrCreate()->AddUpperBound(expr, + IntegerValue(ub)); } PrecedencesPropagator* p = model->GetOrCreate(); diff --git a/ortools/sat/precedences_test.cc b/ortools/sat/precedences_test.cc index 994abc248c..781781d5f2 100644 --- a/ortools/sat/precedences_test.cc +++ b/ortools/sat/precedences_test.cc @@ -45,7 +45,7 @@ using ::testing::UnorderedElementsAre; // TODO(user): move that in a common place. test_utils? #define EXPECT_BOUNDS_EQ(var, lb, ub) \ EXPECT_EQ(integer_trail->LowerBound(var), lb); \ - EXPECT_EQ(integer_trail->UpperBound(var), ub) + EXPECT_EQ(integer_trail->LevelZeroUpperBound(var), ub) // All the tests here uses 10 integer variables initially in [0, 100]. std::vector AddVariables(IntegerTrail* integer_trail) { @@ -68,24 +68,42 @@ TEST(PrecedenceRelationsTest, BasicAPI) { IntegerVariable a(0), b(2), c(4), d(6); PrecedenceRelations precedences(&model); - precedences.Add(a, b, 10); - precedences.Add(d, c, 7); - precedences.Add(b, d, 5); + precedences.AddUpperBound(LinearExpression2::Difference(a, b), -10); + precedences.AddUpperBound(LinearExpression2::Difference(d, c), -7); + precedences.AddUpperBound(LinearExpression2::Difference(b, d), -5); precedences.Build(); - EXPECT_EQ(precedences.GetOffset(a, b), 10); - EXPECT_EQ(precedences.GetOffset(NegationOf(b), NegationOf(a)), 10); - EXPECT_EQ(precedences.GetOffset(a, c), 22); - EXPECT_EQ(precedences.GetOffset(NegationOf(c), NegationOf(a)), 22); - EXPECT_EQ(precedences.GetOffset(a, d), 15); - EXPECT_EQ(precedences.GetOffset(NegationOf(d), NegationOf(a)), 15); - EXPECT_EQ(precedences.GetOffset(d, a), kMinIntegerValue); + EXPECT_EQ( + precedences.LevelZeroUpperBound(LinearExpression2::Difference(a, b)), + -10); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(b), NegationOf(a))), + -10); + EXPECT_EQ( + precedences.LevelZeroUpperBound(LinearExpression2::Difference(a, c)), + -22); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(c), NegationOf(a))), + -22); + EXPECT_EQ( + precedences.LevelZeroUpperBound(LinearExpression2::Difference(a, d)), + -15); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(d), NegationOf(a))), + -15); + EXPECT_EQ( + precedences.LevelZeroUpperBound(LinearExpression2::Difference(d, a)), + kMaxIntegerValue); // Once built, we can update the offsets. // Note however that this would not propagate through the precedence graphs. - precedences.Add(a, b, 15); - EXPECT_EQ(precedences.GetOffset(a, b), 15); - EXPECT_EQ(precedences.GetOffset(NegationOf(b), NegationOf(a)), 15); + precedences.AddUpperBound(LinearExpression2::Difference(a, b), -15); + EXPECT_EQ( + precedences.LevelZeroUpperBound(LinearExpression2::Difference(a, b)), + -15); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(b), NegationOf(a))), + -15); } TEST(PrecedenceRelationsTest, CornerCase1) { @@ -97,16 +115,22 @@ TEST(PrecedenceRelationsTest, CornerCase1) { IntegerVariable a(0), b(2), c(4), d(6); PrecedenceRelations precedences(&model); - precedences.Add(a, b, 10); - precedences.Add(b, c, 7); - precedences.Add(b, d, 5); - precedences.Add(NegationOf(b), a, 5); + precedences.AddUpperBound(LinearExpression2::Difference(a, b), -10); + precedences.AddUpperBound(LinearExpression2::Difference(b, c), -7); + precedences.AddUpperBound(LinearExpression2::Difference(b, d), -5); + precedences.AddUpperBound(LinearExpression2::Difference(NegationOf(b), a), + -5); precedences.Build(); - EXPECT_EQ(precedences.GetOffset(NegationOf(b), a), 5); - EXPECT_EQ(precedences.GetOffset(NegationOf(b), b), 15); - EXPECT_EQ(precedences.GetOffset(NegationOf(b), c), 22); - EXPECT_EQ(precedences.GetOffset(NegationOf(b), d), 20); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(b), a)), + -5); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(b), c)), + -22); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(b), d)), + -20); } TEST(PrecedenceRelationsTest, CornerCase2) { @@ -118,16 +142,34 @@ TEST(PrecedenceRelationsTest, CornerCase2) { IntegerVariable a(0), b(2), c(4), d(6); PrecedenceRelations precedences(&model); - precedences.Add(NegationOf(a), a, 10); - precedences.Add(a, b, 7); - precedences.Add(a, c, 5); - precedences.Add(a, d, 2); + precedences.AddUpperBound(LinearExpression2::Difference(NegationOf(a), a), + -10); + precedences.AddUpperBound(LinearExpression2::Difference(a, b), -7); + precedences.AddUpperBound(LinearExpression2::Difference(a, c), -5); + precedences.AddUpperBound(LinearExpression2::Difference(a, d), -2); + EXPECT_EQ(precedences.LevelZeroUpperBound( + LinearExpression2::Difference(NegationOf(b), NegationOf(a))), + -7); + + precedences.Build(); +} + +TEST(PrecedenceRelationsTest, CoefficientGreaterThanOne) { + Model model; + IntegerTrail* integer_trail = model.GetOrCreate(); + const std::vector vars = AddVariables(integer_trail); + + // Note that odd indices are for the negation. + IntegerVariable a(0), b(2), c(4); + + PrecedenceRelations precedences(&model); + precedences.AddUpperBound(LinearExpression2(a, b, 3, -4), 7); + precedences.AddUpperBound(LinearExpression2(a, c, 2, -3), -5); + precedences.AddUpperBound(LinearExpression2(a, b, 6, -8), 5); + EXPECT_EQ(precedences.LevelZeroUpperBound(LinearExpression2(a, b, 9, -12)), + 6); precedences.Build(); - EXPECT_EQ(precedences.GetOffset(NegationOf(a), a), 10); - EXPECT_EQ(precedences.GetOffset(NegationOf(a), b), 17); - EXPECT_EQ(precedences.GetOffset(NegationOf(a), c), 15); - EXPECT_EQ(precedences.GetOffset(NegationOf(a), d), 12); } TEST(PrecedenceRelationsTest, ConditionalRelations) { @@ -142,20 +184,31 @@ TEST(PrecedenceRelationsTest, ConditionalRelations) { // Note that odd indices are for the negation. IntegerVariable a(0), b(2); PrecedenceRelations precedences(&model); - precedences.PushConditionalRelation({l}, a, b, 15); - precedences.PushConditionalRelation({l}, a, b, 20); + precedences.PushConditionalRelation({l}, LinearExpression2(a, b, 1, 1), 15); + precedences.PushConditionalRelation({l}, LinearExpression2(a, b, 1, 1), 20); // We only keep the best one. - EXPECT_EQ(precedences.GetConditionalOffset(a, NegationOf(b)), -15); - EXPECT_THAT(precedences.GetConditionalEnforcements(a, NegationOf(b)), - ElementsAre(l)); + EXPECT_EQ( + precedences.UpperBound(LinearExpression2::Difference(a, NegationOf(b))), + 15); + std::vector literal_reason; + std::vector integer_reason; + precedences.AddReasonForUpperBoundLowerThan( + LinearExpression2::Difference(a, NegationOf(b)), 15, &literal_reason, + &integer_reason); + EXPECT_THAT(literal_reason, ElementsAre(l.Negated())); // Backtrack works. EXPECT_TRUE(sat_solver->ResetToLevelZero()); - EXPECT_EQ(precedences.GetConditionalOffset(a, NegationOf(b)), - kMinIntegerValue); - EXPECT_THAT(precedences.GetConditionalEnforcements(a, NegationOf(b)), - ElementsAre()); + EXPECT_EQ( + precedences.UpperBound(LinearExpression2::Difference(a, NegationOf(b))), + kMaxIntegerValue); + literal_reason.clear(); + integer_reason.clear(); + precedences.AddReasonForUpperBoundLowerThan( + LinearExpression2::Difference(a, NegationOf(b)), kMaxIntegerValue, + &literal_reason, &integer_reason); + EXPECT_THAT(literal_reason, IsEmpty()); } TEST(PrecedencesPropagatorTest, Empty) { @@ -425,12 +478,18 @@ TEST(PrecedenceRelationsTest, CollectPrecedences) { auto* relations = model.GetOrCreate(); std::vector vars = AddVariables(integer_trail); - relations->Add(vars[0], vars[2], IntegerValue(1)); - relations->Add(vars[0], vars[5], IntegerValue(1)); - relations->Add(vars[1], vars[2], IntegerValue(1)); - relations->Add(vars[2], vars[4], IntegerValue(1)); - relations->Add(vars[3], vars[4], IntegerValue(1)); - relations->Add(vars[4], vars[5], IntegerValue(1)); + relations->AddUpperBound(LinearExpression2::Difference(vars[0], vars[2]), + IntegerValue(-1)); + relations->AddUpperBound(LinearExpression2::Difference(vars[0], vars[5]), + IntegerValue(-1)); + relations->AddUpperBound(LinearExpression2::Difference(vars[1], vars[2]), + IntegerValue(-1)); + relations->AddUpperBound(LinearExpression2::Difference(vars[2], vars[4]), + IntegerValue(-1)); + relations->AddUpperBound(LinearExpression2::Difference(vars[3], vars[4]), + IntegerValue(-1)); + relations->AddUpperBound(LinearExpression2::Difference(vars[4], vars[5]), + IntegerValue(-1)); std::vector p; relations->CollectPrecedences({vars[0], vars[2], vars[3]}, &p); @@ -979,6 +1038,81 @@ TEST(PrecedencesPropagatorTest, BasicFiltering2) { EXPECT_THAT(precedences[1].indices, ElementsAre(0, 1, 2, 3)); } +TEST(BinaryRelationMapsTest, AffineUpperBound) { + Model model; + const IntegerVariable x = model.Add(NewIntegerVariable(0, 10)); + const IntegerVariable y = model.Add(NewIntegerVariable(0, 10)); + const IntegerVariable z = model.Add(NewIntegerVariable(0, 2)); + + // x - y; + LinearExpression2 expr; + expr.vars[0] = x; + expr.vars[1] = y; + expr.coeffs[0] = IntegerValue(1); + expr.coeffs[1] = IntegerValue(-1); + + // Starts with trivial level zero bound. + auto* tested = model.GetOrCreate(); + EXPECT_EQ(tested->UpperBound(expr), IntegerValue(10)); + + // Lets add a relation. + tested->AddRelationBounds(expr, IntegerValue(-5), IntegerValue(5)); + EXPECT_EQ(tested->UpperBound(expr), IntegerValue(5)); + + // Note that we canonicalize with gcd. + expr.coeffs[0] *= 3; + expr.coeffs[1] *= 3; + EXPECT_EQ(tested->UpperBound(expr), IntegerValue(15)); + + // Lets add an affine upper bound to that expression <= 4 * z + 1. + EXPECT_TRUE(tested->AddAffineUpperBound( + expr, AffineExpression(z, IntegerValue(4), IntegerValue(1)))); + EXPECT_EQ(tested->UpperBound(expr), IntegerValue(9)); + + // Lets test the reason, first push a new bound. + auto* search = model.GetOrCreate(); + search->TakeDecision( + Literal(search->GetDecisionLiteral(BooleanOrIntegerLiteral( + IntegerLiteral::LowerOrEqual(z, IntegerValue(1)))))); + + // Because of gcd, even though ub(affine) is now 5, we get 3, + EXPECT_EQ(tested->UpperBound(expr), IntegerValue(3)); + { + std::vector literal_reason; + std::vector integer_reason; + tested->AddReasonForUpperBoundLowerThan(expr, IntegerValue(4), + &literal_reason, &integer_reason); + EXPECT_THAT(literal_reason, ElementsAre()); + EXPECT_THAT(integer_reason, + ElementsAre(IntegerLiteral::LowerOrEqual(z, IntegerValue(1)))); + } + + // If we use a bound not as strong, we get a different reason though. + { + std::vector literal_reason; + std::vector integer_reason; + tested->AddReasonForUpperBoundLowerThan(expr, IntegerValue(9), + &literal_reason, &integer_reason); + EXPECT_THAT(literal_reason, ElementsAre()); + EXPECT_THAT(integer_reason, + ElementsAre(IntegerLiteral::LowerOrEqual(z, IntegerValue(2)))); + } + { + // This is implied by the level zero relation x <= 5 + std::vector literal_reason; + std::vector integer_reason; + tested->AddReasonForUpperBoundLowerThan(expr, IntegerValue(15), + &literal_reason, &integer_reason); + EXPECT_THAT(literal_reason, ElementsAre()); + EXPECT_THAT(integer_reason, ElementsAre()); + } + + // Note that the bound works on the canonicalized expr. + expr.coeffs[0] /= 3; + expr.coeffs[1] /= 3; + EXPECT_EQ(tested->UpperBound(expr), IntegerValue(1)); +} + } // namespace } // namespace sat } // namespace operations_research diff --git a/ortools/sat/sat_inprocessing.cc b/ortools/sat/sat_inprocessing.cc index 4d17f1400c..d25c7897b3 100644 --- a/ortools/sat/sat_inprocessing.cc +++ b/ortools/sat/sat_inprocessing.cc @@ -1636,6 +1636,9 @@ bool BoundedVariableElimination::CrossProduct(BooleanVariable var) { if (new_score_ > score_threshold_) return true; // Perform BVE. + // + // TODO(user): If filter_sat_postsolve_clauses is true, only one of the two + // sets need to be kept for postsolve. if (new_score_ > 0) { if (!ResolveAllClauseContaining(lit)) { diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index c42f5a405e..fb7d23f541 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -24,7 +24,7 @@ option java_multiple_files = true; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 323 +// NEXT TAG: 325 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -411,6 +411,15 @@ message SatParameters { // occurrences of not(x) is not greater than this parameter. optional int32 presolve_bve_threshold = 54 [default = 500]; + // Internal parameter. During BVE, if we eliminate a variable x, by default we + // will push all clauses containing x and all clauses containing not(x) to the + // postsolve. However, it is possible to write the postsolve code so that only + // one such set is needed. The idea is that, if we push the set containing a + // literal l, is to set l to false except if it is needed to satisfy one of + // the clause in the set. This is always beneficial, but for historical + // reason, not all our postsolve algorithm support this. + optional bool filter_sat_postsolve_clauses = 324 [default = false]; + // During presolve, we apply BVE only if this weight times the number of // clauses plus the number of clause literals is not increased. optional int32 presolve_bve_clause_weight = 55 [default = 3]; @@ -937,6 +946,15 @@ message SatParameters { optional int32 maximum_regions_to_split_in_disconnected_no_overlap_2d = 315 [default = 0]; + // When set, this activates a propagator for the no_overlap_2d constraint that + // uses any eventual linear constraints of the model in the form + // `{start interval 1} - {end interval 2} + c*w <= ub` to detect that two + // intervals must overlap in one dimension for some values of `w`. This is + // particularly useful for problems where the distance between two boxes is + // part of the model. + optional bool use_linear3_for_no_overlap_2d_precedences = 323 + [default = true]; + // When set, it activates a few scheduling parameters to improve the lower // bound of scheduling problems. This is only effective with multiple workers // as it modifies the reduced_cost, lb_tree_search, and probing workers. diff --git a/ortools/sat/scheduling_cuts.cc b/ortools/sat/scheduling_cuts.cc index 5fd29afe20..f99d3f0889 100644 --- a/ortools/sat/scheduling_cuts.cc +++ b/ortools/sat/scheduling_cuts.cc @@ -18,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -1121,19 +1120,21 @@ void CtExhaustiveHelper::Init( for (const auto& e2 : events) { if (e2.task_index == e1.task_index) continue; - if (binary_relations->GetPrecedenceStatus(e2.end, e1.start) == + if (binary_relations->GetLevelZeroPrecedenceStatus(e2.end, e1.start) == RelationStatus::IS_TRUE) { predecessors_.AppendToLastVector(e2.task_index); } } } VLOG(2) << "num_tasks:" << max_task_index_ + 1 - << " num_precedences:" << predecessors_.num_entries(); + << " num_precedences:" << predecessors_.num_entries() + << " predecessors size:" << predecessors_.size(); } bool CtExhaustiveHelper::PermutationIsCompatibleWithPrecedences( absl::Span events, absl::Span permutation) { + if (predecessors_.num_entries() == 0) return true; visited_.assign(max_task_index_ + 1, false); for (int i = permutation.size() - 1; i >= 0; --i) { const CompletionTimeEvent& event = events[permutation[i]]; @@ -1159,11 +1160,11 @@ bool ComputeWeightedSumOfEndMinsOfOnePermutationForNoOverlap( IntegerValue end_min_of_previous_task = kMinIntegerValue; for (const int index : permutation) { const CompletionTimeEvent& event = events[index]; - const IntegerValue threshold = + const IntegerValue task_start_min = std::max(event.start_min, end_min_of_previous_task); - if (event.start_max < threshold) return false; // Infeasible. + if (event.start_max < task_start_min) return false; // Infeasible. - end_min_of_previous_task = threshold + event.size_min; + end_min_of_previous_task = task_start_min + event.size_min; sum_of_ends += end_min_of_previous_task; sum_of_weighted_ends += event.energy_min * end_min_of_previous_task; } @@ -1232,14 +1233,13 @@ bool ComputeWeightedSumOfEndMinsOfOnePermutation( // Iterate on the profile to find the step that contains start_min. // Then push until we find a step with enough capacity. - int current = 0; - while (helper.profile_[current + 1].first <= start_min || - helper.profile_[current].second < event.demand_min) { - ++current; + auto profile_it = helper.profile_.begin(); + while ((profile_it + 1)->first <= start_min || + profile_it->second < event.demand_min) { + ++profile_it; } - IntegerValue actual_start = - std::max(start_min, helper.profile_[current].first); + IntegerValue actual_start = std::max(start_min, profile_it->first); const IntegerValue initial_start_min = actual_start; // Propagate precedences. @@ -1255,6 +1255,8 @@ bool ComputeWeightedSumOfEndMinsOfOnePermutation( if (actual_start > initial_start_min) { cut_use_precedences = true; + // Catch up the position on the profile w.r.t. the actual start. + while ((profile_it + 1)->first <= actual_start) ++profile_it; VLOG(3) << "push from " << initial_start_min << " to " << actual_start; } @@ -1276,79 +1278,105 @@ bool ComputeWeightedSumOfEndMinsOfOnePermutation( // Update the profile. helper.new_profile_.clear(); + const IntegerValue demand_min = event.demand_min; + + // Insert the start of the shifted profile. helper.new_profile_.push_back( - {actual_start, helper.profile_[current].second - event.demand_min}); - ++current; + {actual_start, profile_it->second - demand_min}); + ++profile_it; - while (helper.profile_[current].first < actual_end) { + // Copy and modify the part of the profile impacted by the current event. + while (profile_it->first < actual_end) { helper.new_profile_.push_back( - {helper.profile_[current].first, - helper.profile_[current].second - event.demand_min}); - ++current; + {profile_it->first, profile_it->second - demand_min}); + ++profile_it; } - if (helper.profile_[current].first > actual_end) { - helper.new_profile_.push_back( - {actual_end, helper.new_profile_.back().second + event.demand_min}); - } - while (current < helper.profile_.size()) { - helper.new_profile_.push_back(helper.profile_[current]); - ++current; + // Insert a new event in the profile at the end of the task if needed. + if (profile_it->first > actual_end) { + helper.new_profile_.push_back({actual_end, (profile_it - 1)->second}); } + + // Insert the tail of the current profile. + helper.new_profile_.insert(helper.new_profile_.end(), profile_it, + helper.profile_.end()); + helper.profile_.swap(helper.new_profile_); } return true; } +const int kCtExhaustiveTargetSize = 6; +// This correspond to the number of permutations the system will explore when +// fully exploring all possible sizes and all possible permutations for up to 6 +// tasks, without any precedence. +const int kExplorationLimit = 873; // 1! + 2! + 3! + 4! + 5! + 6! + } // namespace -bool ComputeMinSumOfWeightedEndMins( +CompletionTimeExplorationStatus ComputeMinSumOfWeightedEndMins( absl::Span events, IntegerValue capacity_max, double unweighted_threshold, double weighted_threshold, CtExhaustiveHelper& helper, double& min_sum_of_ends, - double& min_sum_of_weighted_ends, bool& cut_use_precedences) { + double& min_sum_of_weighted_ends, bool& cut_use_precedences, + int& exploration_credit) { // Reset the events based sums. min_sum_of_ends = std::numeric_limits::max(); min_sum_of_weighted_ends = std::numeric_limits::max(); + helper.task_to_index_.assign(helper.max_task_index() + 1, -1); + for (int i = 0; i < events.size(); ++i) { + helper.task_to_index_[events[i].task_index] = i; + } + helper.valid_permutation_iterator_.Reset(events.size()); + for (int i = 0; i < events.size(); ++i) { + const int task_i = events[i].task_index; + for (const int task_j : helper.predecessors()[task_i]) { + const int j = helper.task_to_index_[task_j]; + if (j != -1) { + helper.valid_permutation_iterator_.AddArc(j, i); + } + } + } + if (!helper.valid_permutation_iterator_.Init()) { + return CompletionTimeExplorationStatus::NO_VALID_PERMUTATION; + } - // Local stats. - int num_explored = 0; - int num_pruned = 0; - bool aborted = false; - - std::vector permutation(events.size()); - std::iota(permutation.begin(), permutation.end(), 0); + int num_valid_permutations = 0; do { + if (--exploration_credit < 0) break; + IntegerValue sum_of_ends = 0; IntegerValue sum_of_weighted_ends = 0; - if (!helper.PermutationIsCompatibleWithPrecedences(events, permutation)) { - cut_use_precedences = true; - continue; - } if (ComputeWeightedSumOfEndMinsOfOnePermutation( - events, permutation, capacity_max, helper, sum_of_ends, - sum_of_weighted_ends, cut_use_precedences)) { + events, helper.valid_permutation_iterator_.permutation(), + capacity_max, helper, sum_of_ends, sum_of_weighted_ends, + cut_use_precedences)) { min_sum_of_ends = std::min(ToDouble(sum_of_ends), min_sum_of_ends); min_sum_of_weighted_ends = std::min(ToDouble(sum_of_weighted_ends), min_sum_of_weighted_ends); - num_explored++; + ++num_valid_permutations; + if (min_sum_of_ends <= unweighted_threshold && min_sum_of_weighted_ends <= weighted_threshold) { - aborted = true; break; } - } else { - num_pruned++; } - } while (std::next_permutation(permutation.begin(), permutation.end())); - VLOG(3) << "DP: size=" << events.size() << ", explored = " << num_explored - << ", pruned = " << num_pruned << ", aborted = " << aborted - << ", min_sum_of_end_mins = " << min_sum_of_ends - << ", min_sum_of_weighted_end_mins = " << min_sum_of_weighted_ends - << ", unweighted_threshold = " << unweighted_threshold - << ", weighted_threshold = " << weighted_threshold; - return num_explored > 0; + } while (helper.valid_permutation_iterator_.Increase()); + const CompletionTimeExplorationStatus status = + exploration_credit < 0 ? CompletionTimeExplorationStatus::ABORTED + : num_valid_permutations > 0 + ? CompletionTimeExplorationStatus::FINISHED + : CompletionTimeExplorationStatus::NO_VALID_PERMUTATION; + VLOG(2) << "DP: size:" << events.size() + << ", num_valid_permutations:" << num_valid_permutations + << ", min_sum_of_end_mins:" << min_sum_of_ends + << ", min_sum_of_weighted_end_mins:" << min_sum_of_weighted_ends + << ", unweighted_threshold:" << unweighted_threshold + << ", weighted_threshold:" << weighted_threshold + << ", exploration_credit:" << exploration_credit + << ", status:" << status; + return status; } // TODO(user): Improve performance @@ -1359,6 +1387,7 @@ ABSL_MUST_USE_RESULT bool GenerateShortCompletionTimeCutsWithExactBound( IntegerValue capacity_max, CtExhaustiveHelper& helper, Model* model, LinearConstraintManager* manager) { TopNCuts top_n_cuts(5); + // Sort by start min to bucketize by start_min. std::sort( events.begin(), events.end(), @@ -1374,20 +1403,19 @@ ABSL_MUST_USE_RESULT bool GenerateShortCompletionTimeCutsWithExactBound( bool cut_use_precedences = false; // Used for naming the cut. const IntegerValue sequence_start_min = events[start].start_min; - std::vector residual_tasks(events.begin() + start, - events.end()); + helper.residual_events_.assign(events.begin() + start, events.end()); // We look at events that start before sequence_start_min, but are forced // to cross this time point. for (int before = 0; before < start; ++before) { if (events[before].start_min + events[before].size_min > sequence_start_min) { - residual_tasks.push_back(events[before]); // Copy. - residual_tasks.back().lifted = true; + helper.residual_events_.push_back(events[before]); // Copy. + helper.residual_events_.back().lifted = true; } } - std::sort(residual_tasks.begin(), residual_tasks.end(), + std::sort(helper.residual_events_.begin(), helper.residual_events_.end(), [](const CompletionTimeEvent& e1, const CompletionTimeEvent& e2) { return std::tie(e1.lp_end, e1.task_index) < std::tie(e2.lp_end, e2.task_index); @@ -1399,9 +1427,11 @@ ABSL_MUST_USE_RESULT bool GenerateShortCompletionTimeCutsWithExactBound( double sum_of_square_energies = 0; double min_sum_of_ends = std::numeric_limits::max(); double min_sum_of_weighted_ends = std::numeric_limits::max(); + int exploration_limit = kExplorationLimit; + const int kMaxSize = std::min(helper.residual_events_.size(), 12); - for (int i = 0; i < std::min(residual_tasks.size(), 7); ++i) { - const CompletionTimeEvent& event = residual_tasks[i]; + for (int i = 0; i < kMaxSize; ++i) { + const CompletionTimeEvent& event = helper.residual_events_[i]; const double energy = ToDouble(event.energy_min); sum_of_ends_lp += event.lp_end; sum_of_weighted_ends_lp += event.lp_end * energy; @@ -1409,18 +1439,25 @@ ABSL_MUST_USE_RESULT bool GenerateShortCompletionTimeCutsWithExactBound( sum_of_square_energies += energy * energy; // Both cases with 1 or 2 tasks are trivial and independent of the order. - // Also, if capacity is not exceeded, pushing all ends left is a valid LP - // assignment. + // Also, if the sum of demands is less than or equal to the capacity, + // pushing all ends left is a valid LP assignment. And this assignment + // should be propagated by the lp model. if (i <= 1 || sum_of_demands <= capacity_max) continue; - if (!ComputeMinSumOfWeightedEndMins( - absl::MakeSpan(residual_tasks).first(i + 1), capacity_max, + absl::Span tasks_to_explore = + absl::MakeSpan(helper.residual_events_).first(i + 1); + const CompletionTimeExplorationStatus status = + ComputeMinSumOfWeightedEndMins( + tasks_to_explore, capacity_max, /* unweighted_threshold= */ sum_of_ends_lp + kMinCutViolation, /* weighted_threshold= */ sum_of_weighted_ends_lp + kMinCutViolation, helper, min_sum_of_ends, min_sum_of_weighted_ends, - cut_use_precedences)) { + cut_use_precedences, exploration_limit); + if (status == CompletionTimeExplorationStatus::NO_VALID_PERMUTATION) { return false; + } else if (status == CompletionTimeExplorationStatus::ABORTED) { + break; } const double unweigthed_violation = @@ -1435,7 +1472,7 @@ ABSL_MUST_USE_RESULT bool GenerateShortCompletionTimeCutsWithExactBound( LinearConstraintBuilder cut(model, min_sum_of_ends, kMaxIntegerValue); bool is_lifted = false; for (int j = 0; j <= i; ++j) { - const CompletionTimeEvent& event = residual_tasks[j]; + const CompletionTimeEvent& event = helper.residual_events_[j]; is_lifted |= event.lifted; cut.AddTerm(event.end, IntegerValue(1)); } @@ -1452,7 +1489,7 @@ ABSL_MUST_USE_RESULT bool GenerateShortCompletionTimeCutsWithExactBound( kMaxIntegerValue); bool is_lifted = false; for (int j = 0; j <= i; ++j) { - const CompletionTimeEvent& event = residual_tasks[j]; + const CompletionTimeEvent& event = helper.residual_events_[j]; is_lifted |= event.lifted; cut.AddTerm(event.end, event.energy_min); } @@ -1582,8 +1619,7 @@ void AddEventDemandsToCapacitySubsetSum( // ordered by increasing end time in the LP relaxation. void GenerateCompletionTimeCutsWithEnergy( absl::string_view cut_name, std::vector events, - IntegerValue capacity_max, bool skip_low_sizes, Model* model, - LinearConstraintManager* manager) { + IntegerValue capacity_max, Model* model, LinearConstraintManager* manager) { TopNCuts top_n_cuts(5); const VariablesAssignment& assignment = model->GetOrCreate()->Assignment(); @@ -1623,6 +1659,10 @@ void GenerateCompletionTimeCutsWithEnergy( } } + // If we have less than kCtExhaustiveTargetSize tasks, we are already + // covered by the exhaustive cut generator. + if (residual_tasks.size() <= kCtExhaustiveTargetSize) continue; + // Best cut so far for this loop. int best_end = -1; double best_efficacy = 0.01; @@ -1676,9 +1716,8 @@ void GenerateCompletionTimeCutsWithEnergy( AddEventDemandsToCapacitySubsetSum(event, assignment, capacity_max, tmp_possible_demands, dp); - // This is competing with the brute force approach. Skip cases covered - // by the other code. - if (skip_low_sizes && i < 7) continue; + // Ignore cuts covered by the exhaustive cut generator. + if (i < kCtExhaustiveTargetSize) continue; const IntegerValue reachable_capacity = dp.CurrentMax(); @@ -1769,18 +1808,15 @@ CutGenerator CreateNoOverlapCompletionTimeCutGenerator( CtExhaustiveHelper helper; helper.Init(events, model); - const std::string mirror_str = time_is_forward ? "" : "_mirror"; if (!GenerateShortCompletionTimeCutsWithExactBound( - absl::StrCat("NoOverlapCompletionTimeExhaustive", mirror_str), - events, + "NoOverlapCompletionTimeExhaustive", events, /*capacity_max=*/IntegerValue(1), helper, model, manager)) { return false; } GenerateCompletionTimeCutsWithEnergy( - absl::StrCat("NoOverlapCompletionTimeQueyrane", mirror_str), - std::move(events), /*capacity_max=*/IntegerValue(1), - /*skip_low_sizes=*/true, model, manager); + "NoOverlapCompletionTimeQueyrane", std::move(events), + /*capacity_max=*/IntegerValue(1), model, manager); return true; }; if (!generate_cuts(/*time_is_forward=*/true)) return false; @@ -1835,17 +1871,15 @@ CutGenerator CreateCumulativeCompletionTimeCutGenerator( helper.Init(events, model); const IntegerValue capacity_max = integer_trail->UpperBound(capacity); - const std::string mirror_str = time_is_forward ? "" : "_mirror"; if (!GenerateShortCompletionTimeCutsWithExactBound( - absl::StrCat("CumulativeCompletionTimeExhaustive", mirror_str), - events, capacity_max, helper, model, manager)) { + "CumulativeCompletionTimeExhaustive", events, capacity_max, + helper, model, manager)) { return false; } - GenerateCompletionTimeCutsWithEnergy( - absl::StrCat("CumulativeCompletionTimeQueyrane", mirror_str), - std::move(events), capacity_max, - /*skip_low_sizes=*/true, model, manager); + GenerateCompletionTimeCutsWithEnergy("CumulativeCompletionTimeQueyrane", + std::move(events), capacity_max, + model, manager); return true; }; diff --git a/ortools/sat/scheduling_cuts.h b/ortools/sat/scheduling_cuts.h index a1fbd5a61c..920f5a23e6 100644 --- a/ortools/sat/scheduling_cuts.h +++ b/ortools/sat/scheduling_cuts.h @@ -152,6 +152,8 @@ struct CompletionTimeEvent { class CtExhaustiveHelper { public: + CtExhaustiveHelper() = default; + int max_task_index() const { return max_task_index_; } const CompactVectorVector& predecessors() const { return predecessors_; } @@ -159,6 +161,9 @@ class CtExhaustiveHelper { std::vector> profile_; std::vector> new_profile_; std::vector assigned_ends_; + std::vector task_to_index_; + DagTopologicalSortIterator valid_permutation_iterator_; + std::vector residual_events_; // Collect precedences, set max_task_index. // TODO(user): Do some transitive closure. @@ -174,6 +179,27 @@ class CtExhaustiveHelper { std::vector visited_; }; +enum class CompletionTimeExplorationStatus { + FINISHED, + ABORTED, + NO_VALID_PERMUTATION, +}; + +template +void AbslStringify(Sink& sink, const CompletionTimeExplorationStatus& status) { + switch (status) { + case CompletionTimeExplorationStatus::FINISHED: + sink.Append("FINISHED"); + break; + case CompletionTimeExplorationStatus::ABORTED: + sink.Append("ABORTED"); + break; + case CompletionTimeExplorationStatus::NO_VALID_PERMUTATION: + sink.Append("NO_VALID_PERMUTATION"); + break; + } +} + // Computes the minimum sum of the end min and the minimum sum of the end min // weighted by weight of all events. It returns false if no permutation is // valid w.r.t. the range of starts. @@ -182,11 +208,12 @@ class CtExhaustiveHelper { // small, like <= 10. They should also starts in index order. // // Optim: If both sums are proven <= to the corresponding threshold, we abort. -bool ComputeMinSumOfWeightedEndMins( +CompletionTimeExplorationStatus ComputeMinSumOfWeightedEndMins( absl::Span events, IntegerValue capacity_max, double unweighted_threshold, double weighted_threshold, CtExhaustiveHelper& helper, double& min_sum_of_ends, - double& min_sum_of_weighted_ends, bool& cut_use_precedences); + double& min_sum_of_weighted_ends, bool& cut_use_precedences, + int& exploration_credit); } // namespace sat } // namespace operations_research diff --git a/ortools/sat/scheduling_cuts_test.cc b/ortools/sat/scheduling_cuts_test.cc index f434e296de..543bafd019 100644 --- a/ortools/sat/scheduling_cuts_test.cc +++ b/ortools/sat/scheduling_cuts_test.cc @@ -411,9 +411,12 @@ TEST(ComputeMinSumOfEndMinsTest, CombinationOf3) { CtExhaustiveHelper ct_helper; ct_helper.Init(events, &model); bool cut_use_precedences = false; - ASSERT_TRUE(ComputeMinSumOfWeightedEndMins( - events, two, 0.01, 0.01, ct_helper, min_sum_of_end_mins, - min_sum_of_weighted_end_mins, cut_use_precedences)); + int exploration_credit = 1000; + ASSERT_EQ(ComputeMinSumOfWeightedEndMins( + events, two, 0.01, 0.01, ct_helper, min_sum_of_end_mins, + min_sum_of_weighted_end_mins, cut_use_precedences, + exploration_credit), + CompletionTimeExplorationStatus::FINISHED); EXPECT_EQ(min_sum_of_end_mins, 17); EXPECT_EQ(min_sum_of_weighted_end_mins, 86); EXPECT_FALSE(cut_use_precedences); @@ -460,13 +463,67 @@ TEST(ComputeMinSumOfEndMinsTest, CombinationOf3ConstraintStart) { CtExhaustiveHelper ct_helper; ct_helper.Init(events, &model); bool cut_use_precedences = false; - ASSERT_TRUE(ComputeMinSumOfWeightedEndMins( - events, two, 0.01, 0.01, ct_helper, min_sum_of_end_mins, - min_sum_of_weighted_end_mins, cut_use_precedences)); + int exploration_credit = 1000; + + ASSERT_EQ(ComputeMinSumOfWeightedEndMins( + events, two, 0.01, 0.01, ct_helper, min_sum_of_end_mins, + min_sum_of_weighted_end_mins, cut_use_precedences, + exploration_credit), + CompletionTimeExplorationStatus::FINISHED); EXPECT_EQ(min_sum_of_end_mins, 18); EXPECT_EQ(min_sum_of_weighted_end_mins, 86); } +TEST(ComputeMinSumOfEndMinsTest, Abort) { + Model model; + auto* intervals_repository = model.GetOrCreate(); + + IntegerValue one(1); + IntegerValue two(2); + + const IntegerVariable start1 = model.Add(NewIntegerVariable(0, 3)); + const IntegerValue size1(3); + const IntervalVariable i1 = intervals_repository->CreateInterval( + start1, AffineExpression(start1, one, size1), size1, kNoLiteralIndex, + /*add_linear_relation=*/false); + + const IntegerVariable start2 = model.Add(NewIntegerVariable(0, 10)); + const IntegerValue size2(4); + const IntervalVariable i2 = intervals_repository->CreateInterval( + start2, AffineExpression(start2, one, size2), size2, kNoLiteralIndex, + /*add_linear_relation=*/false); + + const IntegerVariable start3 = model.Add(NewIntegerVariable(0, 10)); + const IntegerValue size3(5); + const IntervalVariable i3 = intervals_repository->CreateInterval( + start3, AffineExpression(start3, one, size3), size3, kNoLiteralIndex, + /*add_linear_relation=*/false); + + SchedulingConstraintHelper* helper = + model.GetOrCreate()->GetOrCreateHelper({i1, i2, i3}); + SchedulingDemandHelper* demands_helper = + new SchedulingDemandHelper({two, one, one}, helper, &model); + model.TakeOwnership(demands_helper); + + CompletionTimeEvent e1(0, helper, demands_helper); + CompletionTimeEvent e2(1, helper, demands_helper); + CompletionTimeEvent e3(2, helper, demands_helper); + const std::vector events = {e1, e2, e3}; + + double min_sum_of_end_mins = 0; + double min_sum_of_weighted_end_mins = 0; + CtExhaustiveHelper ct_helper; + ct_helper.Init(events, &model); + bool cut_use_precedences = false; + int exploration_credit = 2; + + ASSERT_EQ(ComputeMinSumOfWeightedEndMins( + events, two, 0.01, 0.01, ct_helper, min_sum_of_end_mins, + min_sum_of_weighted_end_mins, cut_use_precedences, + exploration_credit), + CompletionTimeExplorationStatus::ABORTED); +} + TEST(ComputeMinSumOfEndMinsTest, Infeasible) { Model model; auto* intervals_repository = model.GetOrCreate(); @@ -508,9 +565,12 @@ TEST(ComputeMinSumOfEndMinsTest, Infeasible) { CtExhaustiveHelper ct_helper; ct_helper.Init(events, &model); bool cut_use_precedences = false; - ASSERT_FALSE(ComputeMinSumOfWeightedEndMins( - events, two, 0.01, 0.01, ct_helper, min_sum_of_end_mins, - min_sum_of_weighted_end_mins, cut_use_precedences)); + int exploration_credit = 1000; + ASSERT_EQ(ComputeMinSumOfWeightedEndMins( + events, two, 0.01, 0.01, ct_helper, min_sum_of_end_mins, + min_sum_of_weighted_end_mins, cut_use_precedences, + exploration_credit), + CompletionTimeExplorationStatus::NO_VALID_PERMUTATION); } double ExactMakespan(absl::Span sizes, std::vector& demands, @@ -570,15 +630,18 @@ double ExactMakespanBruteForce(absl::Span sizes, CtExhaustiveHelper ct_helper; ct_helper.Init(events, &model); bool cut_use_precedences = false; - EXPECT_TRUE(ComputeMinSumOfWeightedEndMins( - events, capacity, 0.01, 0.01, ct_helper, min_sum_of_end_mins, - min_sum_of_weighted_end_mins, cut_use_precedences)); + int exploration_credit = 10000; + EXPECT_EQ(ComputeMinSumOfWeightedEndMins( + events, capacity, 0.01, 0.01, ct_helper, min_sum_of_end_mins, + min_sum_of_weighted_end_mins, cut_use_precedences, + exploration_credit), + CompletionTimeExplorationStatus::FINISHED); return min_sum_of_end_mins; } TEST(ComputeMinSumOfEndMinsTest, RandomCases) { absl::BitGen random; - const int kNumTests = DEBUG_MODE ? 100 : 1000; + const int kNumTests = DEBUG_MODE ? 50 : 500; const int kNumTasks = 7; for (int loop = 0; loop < kNumTests; ++loop) { const int capacity = absl::Uniform(random, 10, 30); diff --git a/ortools/sat/scheduling_helpers.cc b/ortools/sat/scheduling_helpers.cc index bc596380ec..69d0bad256 100644 --- a/ortools/sat/scheduling_helpers.cc +++ b/ortools/sat/scheduling_helpers.cc @@ -346,26 +346,20 @@ IntegerValue SchedulingConstraintHelper::GetCurrentMinDistanceBetweenTasks( int a, int b, bool add_reason_if_after) { const AffineExpression before = ends_[a]; const AffineExpression after = starts_[b]; - if (before.var == kNoIntegerVariable || before.coeff != 1 || - after.var == kNoIntegerVariable || after.coeff != 1) { - return kMinIntegerValue; - } + LinearExpression2 expr(before.var, after.var, before.coeff, -after.coeff); - // We take the max of the level zero offset and the one coming from a - // conditional precedence at true. - const IntegerValue conditional_offset = - precedence_relations_->GetConditionalOffset(before.var, after.var); - const IntegerValue known = integer_trail_->LevelZeroLowerBound(after.var) - - integer_trail_->LevelZeroUpperBound(before.var); - const IntegerValue offset = std::max(conditional_offset, known); + // We take the min of the level zero (end_a - start_b) and the one coming from + // a conditional precedence at true. + const IntegerValue conditional_ub = precedence_relations_->UpperBound(expr); + const IntegerValue level_zero_ub = integer_trail_->LevelZeroUpperBound(expr); + const IntegerValue expr_ub = std::min(conditional_ub, level_zero_ub); const IntegerValue needed_offset = before.constant - after.constant; - const IntegerValue distance = offset - needed_offset; - if (add_reason_if_after && distance >= 0 && known < conditional_offset) { - for (const Literal l : precedence_relations_->GetConditionalEnforcements( - before.var, after.var)) { - literal_reason_.push_back(l.Negated()); - } + const IntegerValue ub_of_end_minus_start = expr_ub + needed_offset; + const IntegerValue distance = -ub_of_end_minus_start; + if (add_reason_if_after && distance >= 0 && level_zero_ub > conditional_ub) { + precedence_relations_->AddReasonForUpperBoundLowerThan( + expr, conditional_ub, MutableLiteralReason(), MutableIntegerReason()); } return distance; } @@ -394,7 +388,9 @@ bool SchedulingConstraintHelper::PropagatePrecedence(int a, int b) { } } const IntegerValue offset = before.constant - after.constant; - if (precedence_relations_->Add(before.var, after.var, offset)) { + const LinearExpression2 expr = + LinearExpression2::Difference(before.var, after.var); + if (precedence_relations_->AddUpperBound(expr, -offset)) { VLOG(2) << "new relation " << TaskDebugString(a) << " <= " << TaskDebugString(b); if (before.var == NegationOf(after.var)) { diff --git a/ortools/sat/simplification.cc b/ortools/sat/simplification.cc index 12497a1e91..590b2f1be1 100644 --- a/ortools/sat/simplification.cc +++ b/ortools/sat/simplification.cc @@ -231,6 +231,18 @@ void SatPresolver::SetNumVariables(int num_variables) { } } +void SatPresolver::RebuildLiteralToClauses() { + const int size = literal_to_clauses_.size(); + literal_to_clauses_.clear(); + literal_to_clauses_.resize(size); + for (ClauseIndex ci(0); ci < clauses_.size(); ++ci) { + for (const Literal lit : clauses_[ci]) { + literal_to_clauses_[lit].push_back(ci); + } + } + num_deleted_literals_since_last_cleanup_ = 0; +} + void SatPresolver::AddClauseInternal(std::vector* clause) { if (drat_proof_handler_ != nullptr) drat_proof_handler_->AddClause(*clause); @@ -365,6 +377,10 @@ bool SatPresolver::Presolve(const std::vector& can_be_removed) { if (!can_be_removed[var.value()]) continue; if (CrossProduct(Literal(var, true))) { if (!ProcessAllClauses()) return false; + + if (num_deleted_literals_since_last_cleanup_ > 1e7) { + RebuildLiteralToClauses(); + } } if (time_limit_ != nullptr && time_limit_->LimitReached()) return true; if (num_inspected_signatures_ + num_inspected_literals_ > 1e9) return true; @@ -702,9 +718,18 @@ bool SatPresolver::ProcessClauseToSimplifyOthers(ClauseIndex clause_index) { return true; } -void SatPresolver::RemoveAndRegisterForPostsolveAllClauseContaining(Literal x) { - for (ClauseIndex i : literal_to_clauses_[x]) { - if (!clauses_[i].empty()) RemoveAndRegisterForPostsolve(i, x); +void SatPresolver::RemoveAllClauseContaining(Literal x, + bool register_for_postsolve) { + if (register_for_postsolve) { + for (ClauseIndex i : literal_to_clauses_[x]) { + if (!clauses_[i].empty()) { + RemoveAndRegisterForPostsolve(i, x); + } + } + } else { + for (ClauseIndex i : literal_to_clauses_[x]) { + if (!clauses_[i].empty()) Remove(i); + } } gtl::STLClearObject(&literal_to_clauses_[x]); literal_to_clause_sizes_[x] = 0; @@ -725,18 +750,24 @@ bool SatPresolver::CrossProduct(Literal x) { } // Compute the threshold under which we don't remove x.Variable(). - int threshold = 0; + int num_clauses = 0; + int64_t sum_for_x = 0; + int64_t sum_for_not_x = 0; const int clause_weight = parameters_.presolve_bve_clause_weight(); for (ClauseIndex i : literal_to_clauses_[x]) { if (!clauses_[i].empty()) { - threshold += clause_weight + clauses_[i].size(); + ++num_clauses; + sum_for_x += clauses_[i].size(); } } for (ClauseIndex i : literal_to_clauses_[x.NegatedIndex()]) { if (!clauses_[i].empty()) { - threshold += clause_weight + clauses_[i].size(); + ++num_clauses; + sum_for_not_x += clauses_[i].size(); } } + const int64_t threshold = + clause_weight * num_clauses + sum_for_x + sum_for_not_x; // For the BCE, we prefer s2 to be small. if (s1 < s2) x = x.Negated(); @@ -792,8 +823,17 @@ bool SatPresolver::CrossProduct(Literal x) { // // TODO(user): We could only update the priority queue once for each variable // instead of doing it many times. - RemoveAndRegisterForPostsolveAllClauseContaining(x); - RemoveAndRegisterForPostsolveAllClauseContaining(x.Negated()); + bool push_x_for_postsolve = true; + bool push_not_x_for_postsolve = true; + if (parameters_.filter_sat_postsolve_clauses()) { + if (sum_for_x <= sum_for_not_x) { + push_not_x_for_postsolve = false; + } else { + push_x_for_postsolve = false; + } + } + RemoveAllClauseContaining(x, push_x_for_postsolve); + RemoveAllClauseContaining(x.Negated(), push_not_x_for_postsolve); // TODO(user): At this point x.Variable() is added back to the priority queue. // Avoid doing that. @@ -801,8 +841,9 @@ bool SatPresolver::CrossProduct(Literal x) { } void SatPresolver::Remove(ClauseIndex ci) { + num_deleted_literals_since_last_cleanup_ += clauses_[ci].size(); signatures_[ci] = 0; - for (Literal e : clauses_[ci]) { + for (const Literal e : clauses_[ci]) { literal_to_clause_sizes_[e]--; UpdatePriorityQueue(e.Variable()); UpdateBvaPriorityQueue(Literal(e.Variable(), true).Index()); diff --git a/ortools/sat/simplification.h b/ortools/sat/simplification.h index 75627dd2f1..1f29de687d 100644 --- a/ortools/sat/simplification.h +++ b/ortools/sat/simplification.h @@ -244,10 +244,14 @@ class SatPresolver { // after this call. void AddClauseInternal(std::vector* clause); + // Since we only cleanup the list lazily, literal_to_clauses_ memory usage + // can get out of hand, we clean it up periodically. + void RebuildLiteralToClauses(); + // Clause removal function. void Remove(ClauseIndex ci); void RemoveAndRegisterForPostsolve(ClauseIndex ci, Literal x); - void RemoveAndRegisterForPostsolveAllClauseContaining(Literal x); + void RemoveAllClauseContaining(Literal x, bool register_for_postsolve); // Call ProcessClauseToSimplifyOthers() on all the clauses in // clause_to_process_ and empty the list afterwards. Note that while some @@ -355,6 +359,10 @@ class SatPresolver { // Occurrence list. For each literal, contains the ClauseIndex of the clause // that contains it (ordered by clause index). + // + // This is cleaned up lazily, or when num_deleted_literals_since_last_cleanup_ + // becomes big. + int64_t num_deleted_literals_since_last_cleanup_ = 0; util_intops::StrongVector> literal_to_clauses_; diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index e3cd1cd66d..a4272ca360 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -164,7 +164,7 @@ void SharedResponseManager::LogMessage(absl::string_view prefix, } void SharedResponseManager::LogMessageWithThrottling( - const std::string& prefix, const std::string& message) { + absl::string_view prefix, absl::string_view message) { absl::MutexLock mutex_lock(&mutex_); int id; diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index fe6a22fdf8..c6eadff080 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -399,8 +399,8 @@ class SharedResponseManager { // Wrapper around our SolverLogger, but protected by mutex. void LogMessage(absl::string_view prefix, absl::string_view message); - void LogMessageWithThrottling(const std::string& prefix, - const std::string& message); + void LogMessageWithThrottling(absl::string_view prefix, + absl::string_view message); bool LoggingIsEnabled() const; void AppendResponseToBeMerged(const CpSolverResponse& response); diff --git a/ortools/sat/util.h b/ortools/sat/util.h index ea18680bc4..4c2a26b9f3 100644 --- a/ortools/sat/util.h +++ b/ortools/sat/util.h @@ -26,6 +26,7 @@ #include #include +#include "absl/base/attributes.h" #include "absl/container/btree_set.h" #include "absl/log/check.h" #include "absl/log/log_streamer.h" @@ -310,7 +311,7 @@ bool LinearInequalityCanBeReducedWithClosestMultiple( // The model "singleton" random engine used in the solver. // // In test, we usually set use_absl_random() so that the sequence is changed at -// each invocation. This way, clients do not relly on the wrong assumption that +// each invocation. This way, clients do not really on the wrong assumption that // a particular optimal solution will be returned if they are many equivalent // ones. class ModelRandomGenerator : public absl::BitGenRef { @@ -979,6 +980,196 @@ inline void CompactVectorVector::ResetFromTranspose( starts_[0] = 0; } +// A class to generate all possible topological sorting of a dag. +// +// If the graph has no edges, it will generate all possible permutations. +// +// If the graph has edges, it will generate all possible permutations of the +// dag that are a topological sorting of the graph. +// +// The class maintains 5 fields: +// - graph_: a vector of vectors, where graph_[i] contains the list of +// elements that are adjacent to element i. +// - size_: the size of the graph. +// +// - missing_parent_numbers_: a vector of integers, where +// missing_parent_numbers_[i] is the number of parents of element i that are +// not yet in permutation_. Before Init() is called, no element is yet in +// permutation_ so that it is the number of parents of i. After Init(), and +// before Increase() returns true, it is always 0 (except during the +// execution of Increase(), see below). +// +// - permutation_: a vector of integers, that after Init() is called, and +// before Increase() returns false, it is a topological sorting of the graph +// (except during the execution of Increase()). +// - element_original_position_: a vector of integers, where +// element_original_position_[i] is the original position of element i in the +// permutation_. See the algorithm below for more details. + +class DagTopologicalSortIterator { + public: + // Graph maps indices to their children. Any children must exist. + DagTopologicalSortIterator() : size_(0) {} + explicit DagTopologicalSortIterator(int size) : size_(size) { Reset(size); } + + void Reset(int size) { + size_ = size; + graph_.assign(size, {}); + missing_parent_numbers_.assign(size, 0); + permutation_.clear(); + element_original_position_.assign(size, 0); + } + + // Must be called before Init(). + void AddArc(int from, int to) { + DCHECK_GE(from, 0); + DCHECK_LT(from, size_); + DCHECK_GE(to, 0); + DCHECK_LT(to, size_); + graph_[from].push_back(to); + missing_parent_numbers_[to]++; + } + + // To describe the algorithm in Increase() and Init(), we consider the + // following invariant, called Invariant(pos) for a position pos in [0, + // size_): + // 1. permutations_[0], ..., permutations_[pos] form a prefix of a + // topological ordering of the graph; + // 2. permutations_[pos + 1], ..., permutations_.back() are all other + // elements that have all their parents in permutations_[0], ..., + // permutations_[pos], ordered lexicographically by the index of their + // last parent in permutations_[0], ... permutations_[pos] and then by + // their index in the graph; + // 3. missing_parent_numbers_[i] is the number of parents of element i that + // are not in {permutations_[0], ..., permutations_[pos]}. + // 4. element_original_position_[i] is the original position of element i of + // the permutation following the order described in 2. In particular, + // element_original_position_[i] = i for i > pos. + // Set and Unset maintain these invariants. + + // Precondition: Invariant(size_ - 1) holds. + // Postcondition: Invariant(size_ - 1) holds if Increase() returns true. + // If Increase() returns false, all topological orderings of the graph have + // been generated and the state of permutation_ is not specified.. + bool Increase() { + Unset(size_ - 1); + for (int pos = size_ - 2; pos >= 0; --pos) { + // Invariant(pos) holds. + // Increasing logic: once permutation_[pos] has been put back to its + // original position by Unset(pos), elements permutations_[pos], ..., + // permutations_.back() are in their original ordering, in particular in + // the same order as last time the iteration on permutation_[pos] + // occurred (according to Invariant(pos).2, these are exactly the elements + // that have to be tried at pos). + // All possibilities in permutations_[pos], ..., + // permutations_[element_original_position_[pos]] have been run through. + // The next to test is permutations_[element_original_position_[pos] + 1]. + const int k = element_original_position_[pos] + 1; + Unset(pos); + // Invariant(pos - 1) holds. + + // No more elements to iterate on at position pos. Go backwards one + // position to increase that one. + if (k == permutation_.size()) continue; + Set(pos, k); + // Invariant(pos) holds. + for (++pos; pos < size_; ++pos) { + // Invariant(pos - 1) holds. + // According to Invariant(pos - 1).2, if pos >= permutation_.size(), + // there are no more elements we can add to the permutation which means + // that we detected a cycle. It would be a bug as we would have detected + // it in Init(). + CHECK_LT(pos, permutation_.size()) << "Cycle detected"; + // According to Invariant(pos - 1).2, elements that can be used at pos + // are permutations_[pos], ..., permutations_.back(). Starts the + // iteration at permutations_[pos]. + Set(pos, pos); + // Invariant(pos) holds. + } + // Invariant(size_ - 1) holds. + return true; + } + return false; + } + + // Must be called before Increase(). + ABSL_MUST_USE_RESULT bool Init() { + for (int i = 0; i < size_; ++i) { + if (missing_parent_numbers_[i] == 0) { + permutation_.push_back(i); + } + } + for (int pos = 0; pos < size_; ++pos) { + // Invariant(pos - 1) holds. + // According to Invariant(pos - 1).2, if pos >= permutation_.size(), + // there are no more elements we can add to the permutation. + if (pos >= permutation_.size()) return false; + // According to Invariant(pos - 1).2, elements that can be used at pos + // are permutations_[pos], ..., permutations_.back(). Starts the + // iteration at permutations_[pos]. + Set(pos, pos); + // Invariant(pos) holds. + } + // Invariant(pos - 1) hold. We have a permutation. + return true; + } + + const std::vector& permutation() const { return permutation_; } + + private: + // Graph maps indices to their children. Children must be in [0, size_). + std::vector> graph_; + // Number of elements in graph_. + int size_; + // For each element in graph_, the number of parents it has that are not yet + // in permutation_. In particular, it is always 0 when Init has been called + // and when Increase is not in progress (and has not yet returned false). + std::vector missing_parent_numbers_; + // The current permutation. It is ensured to be a topological sorting of the + // graph once Init has been called and Increase has not yet returned false. + std::vector permutation_; + // Keeps track of the original position of the element in permutation_[i]. See + // the comment above the class for the detailed algorithm. + std::vector element_original_position_; + + // Unset the element at pos. + // + // - Precondition: Invariant(pos) holds. + // - Postcondition: Invariant(pos - 1) holds. + void Unset(int pos) { + const int n = permutation_[pos]; + // Before the loop: Invariant(pos).2 and Invariant(pos).3 hold. + // After the swap below: Invariant(pos - 1).2 and Invariant(pos - 1).3 hold. + for (const int c : graph_[n]) { + if (missing_parent_numbers_[c] == 0) permutation_.pop_back(); + ++missing_parent_numbers_[c]; + } + std::swap(permutation_[element_original_position_[pos]], permutation_[pos]); + // Invariant(pos).4 -> Invariant(pos - 1).4. + element_original_position_[pos] = pos; + } + + // Set the element at pos to the element at k. + // + // - Precondition: Invariant(pos - 1) holds and k in [pos, + // permutation_.size()). + // - Postcondition: Invariant(pos) holds and permutation_[pos] has been + // swapped with permutation_[k]. + void Set(int pos, int k) { + int n = permutation_[k]; + // Before the loop: Invariant(pos - 1).2 and Invariant(pos - 1).3 hold. + // After the loop: Invariant(pos).2 and Invariant(pos).3 hold. + for (int c : graph_[n]) { + --missing_parent_numbers_[c]; + if (missing_parent_numbers_[c] == 0) permutation_.push_back(c); + } + // Invariant(pos - 1).1 -> Invariant(pos).1. + std::swap(permutation_[k], permutation_[pos]); + // Invariant(pos - 1).4 -> Invariant(pos).4. + element_original_position_[pos] = k; + } +}; + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/util_test.cc b/ortools/sat/util_test.cc index 17a06e893c..1e06fbfb2d 100644 --- a/ortools/sat/util_test.cc +++ b/ortools/sat/util_test.cc @@ -23,9 +23,11 @@ #include #include #include +#include #include #include "absl/container/btree_set.h" +#include "absl/container/flat_hash_set.h" #include "absl/log/check.h" #include "absl/numeric/int128.h" #include "absl/random/random.h" @@ -150,7 +152,7 @@ TEST(CompactVectorVectorTest, ResetFromTranspose) { EXPECT_THAT(transpose[4], ElementsAre(2)); EXPECT_THAT(transpose[5], ElementsAre(0)); - // Note that retransposing sorts ! + // Note that re-transposing sorts ! CompactVectorVector second_transpose; second_transpose.ResetFromTranspose(transpose); @@ -1065,6 +1067,105 @@ TEST(MaxBoundedSubsetSumExactTest, CornerCase) { EXPECT_EQ(helper.MaxSubsetSum({0, 5, 6}, 4), 0); } +TEST(DagTopologicalSortIteratorTest, GenerateValidPermutations) { + DagTopologicalSortIterator dag_iterator(6); + dag_iterator.AddArc(5, 2); + dag_iterator.AddArc(5, 0); + dag_iterator.AddArc(4, 0); + dag_iterator.AddArc(4, 1); + dag_iterator.AddArc(2, 3); + dag_iterator.AddArc(3, 1); + EXPECT_TRUE(dag_iterator.Init()); + int count = 0; + do { + ++count; + } while (dag_iterator.Increase()); + EXPECT_EQ(count, 13); +} + +TEST(DagTopologicalSortIteratorTest, GenerateAllPermutations) { + DagTopologicalSortIterator dag_iterator(6); + EXPECT_TRUE(dag_iterator.Init()); + int count = 0; + do { + ++count; + } while (dag_iterator.Increase()); + EXPECT_EQ(count, 720); +} + +TEST(DagTopologicalSortIteratorTest, OnePrecedence) { + DagTopologicalSortIterator dag_iterator(6); + dag_iterator.AddArc(5, 2); + EXPECT_TRUE(dag_iterator.Init()); + int count = 0; + do { + ++count; + } while (dag_iterator.Increase()); + EXPECT_EQ(count, 360); +} + +TEST(DagTopologicalSortIteratorTest, ReversePrecedence) { + DagTopologicalSortIterator dag_iterator(6); + dag_iterator.AddArc(2, 5); + EXPECT_TRUE(dag_iterator.Init()); + int count = 0; + do { + ++count; + } while (dag_iterator.Increase()); + EXPECT_EQ(count, 360); +} + +TEST(DagTopologicalSortIteratorTest, RandomTest) { + absl::BitGen random; + for (int i = 0; i < 5000; ++i) { + DagTopologicalSortIterator dag_iterator(6); + + const int num_arcs = absl::Uniform(random, 1, 10); + absl::flat_hash_set> arcs; + while (arcs.size() < num_arcs) { + const int from = absl::Uniform(random, 0, 6); + int to = absl::Uniform(random, 0, 5); + if (from == to) ++to; + if (arcs.insert({from, to}).second) { + dag_iterator.AddArc(from, to); + } + } + + absl::flat_hash_set> iterator_solutions; + int count_iterator = 0; + if (dag_iterator.Init()) { + do { + ++count_iterator; + iterator_solutions.insert(dag_iterator.permutation()); + } while (dag_iterator.Increase()); + } + + std::vector permutation = {0, 1, 2, 3, 4, 5}; + absl::flat_hash_set> permutation_solutions; + int count_permutation = 0; + do { + bool ok = true; + for (int i = 1; i < permutation.size(); ++i) { + if (!ok) break; + const int after = permutation[i]; + for (int j = 0; j < i; ++j) { + const int before = permutation[j]; + if (arcs.contains({after, before})) { + ok = false; + break; + } + } + } + if (ok) { + ++count_permutation; + permutation_solutions.insert(permutation); + } + } while (std::next_permutation(permutation.begin(), permutation.end())); + EXPECT_EQ(count_permutation, count_iterator); + EXPECT_EQ(iterator_solutions, permutation_solutions); + } +} + } // namespace } // namespace sat } // namespace operations_research