diff --git a/ortools/lp_data/sparse_vector.h b/ortools/lp_data/sparse_vector.h index 5dda12520e..571cccf9df 100644 --- a/ortools/lp_data/sparse_vector.h +++ b/ortools/lp_data/sparse_vector.h @@ -503,8 +503,8 @@ SparseVector::SparseVector(const SparseVector& other) { } template -SparseVector& SparseVector:: -operator=(const SparseVector& other) { +SparseVector& +SparseVector::operator=(const SparseVector& other) { PopulateFromSparseVector(other); return *this; } diff --git a/ortools/sat/all_different.cc b/ortools/sat/all_different.cc index 7c642837d6..e325f496c1 100644 --- a/ortools/sat/all_different.cc +++ b/ortools/sat/all_different.cc @@ -56,7 +56,8 @@ std::function AllDifferentBinary( std::function AllDifferentOnBounds( const std::vector& vars) { return [=](Model* model) { - AllDifferentBoundsPropagator* constraint = new AllDifferentBoundsPropagator( + if (vars.empty()) return; + auto* constraint = new AllDifferentBoundsPropagator( vars, model->GetOrCreate()); constraint->RegisterWith(model->GetOrCreate()); model->TakeOwnership(constraint); @@ -407,14 +408,22 @@ bool AllDifferentConstraint::Propagate() { AllDifferentBoundsPropagator::AllDifferentBoundsPropagator( const std::vector& vars, IntegerTrail* integer_trail) - : vars_(vars), integer_trail_(integer_trail), num_calls_(0) { + : integer_trail_(integer_trail) { + CHECK(!vars.empty()); + + // We need +2 for sentinels. + const int capacity = vars.size() + 2; + index_to_start_index_.resize(capacity); + index_to_end_index_.resize(capacity); + index_to_var_.resize(capacity, kNoIntegerVariable); + for (int i = 0; i < vars.size(); ++i) { - negated_vars_.push_back(NegationOf(vars_[i])); + vars_.push_back({vars[i]}); + negated_vars_.push_back({NegationOf(vars[i])}); } } bool AllDifferentBoundsPropagator::Propagate() { - if (vars_.empty()) return true; if (!PropagateLowerBounds()) return false; // Note that it is not required to swap back vars_ and negated_vars_. @@ -425,77 +434,183 @@ bool AllDifferentBoundsPropagator::Propagate() { return result; } -// TODO(user): we could gain by pushing all the new bound at the end, so that -// we just have to sort to_insert_ once. void AllDifferentBoundsPropagator::FillHallReason(IntegerValue hall_lb, IntegerValue hall_ub) { - for (auto entry : to_insert_) { - value_to_variable_[entry.first] = entry.second; - } - to_insert_.clear(); integer_reason_.clear(); - for (int64 v = hall_lb.value(); v <= hall_ub; ++v) { - const IntegerVariable var = gtl::FindOrDie(value_to_variable_, v); + const int limit = GetIndex(hall_ub); + for (int i = GetIndex(hall_lb); i <= limit; ++i) { + const IntegerVariable var = index_to_var_[i]; integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, hall_lb)); integer_reason_.push_back(IntegerLiteral::LowerOrEqual(var, hall_ub)); } } -bool AllDifferentBoundsPropagator::PropagateLowerBounds() { - ++num_calls_; - critical_intervals_.clear(); - hall_starts_.clear(); - hall_ends_.clear(); - - to_insert_.clear(); - if (num_calls_ % 20 == 0) { - // We don't really need to clear this, but we do from time to time to - // save memory (in case the variable domains are huge). This optimization - // helps a bit. - value_to_variable_.clear(); +int AllDifferentBoundsPropagator::FindStartIndexAndCompressPath(int index) { + // First, walk the pointer until we find one pointing to itself. + int start_index = index; + while (true) { + const int next = index_to_start_index_[start_index]; + if (start_index == next) break; + start_index = next; } - // Loop over the variables by increasing ub. - IncrementalSort( - vars_.begin(), vars_.end(), [this](IntegerVariable a, IntegerVariable b) { - return integer_trail_->UpperBound(a) < integer_trail_->UpperBound(b); - }); - for (const IntegerVariable var : vars_) { - const IntegerValue lb = integer_trail_->LowerBound(var); + // Second, redo the same thing and make everyone point to the representative. + while (true) { + const int next = index_to_start_index_[index]; + if (start_index == next) break; + index_to_start_index_[index] = start_index; + index = next; + } + return start_index; +} - // Check if lb is in an Hall interval, and push it if this is the case. - const int hall_index = - std::lower_bound(hall_ends_.begin(), hall_ends_.end(), lb) - - hall_ends_.begin(); - if (hall_index < hall_ends_.size() && hall_starts_[hall_index] <= lb) { - const IntegerValue hs = hall_starts_[hall_index]; - const IntegerValue he = hall_ends_[hall_index]; - FillHallReason(hs, he); - integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, hs)); - if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(var, he + 1), - /*literal_reason=*/{}, integer_reason_)) { +bool AllDifferentBoundsPropagator::PropagateLowerBounds() { + // Start by filling the cached bounds and sorting by increasing lb. + for (VarValue& entry : vars_) { + entry.lb = integer_trail_->LowerBound(entry.var); + entry.ub = integer_trail_->UpperBound(entry.var); + } + IncrementalSort(vars_.begin(), vars_.end(), + [](VarValue a, VarValue b) { return a.lb < b.lb; }); + + // We will split the variable in vars sorted by lb in contiguous subset with + // index of the form [start, start + num_in_window). + int start = 0; + int num_in_window = 1; + + // Minimum lower bound in the current window. + IntegerValue min_lb = vars_.front().lb; + + const int size = vars_.size(); + for (int i = 1; i < size; ++i) { + const IntegerValue lb = vars_[i].lb; + + // If the lower bounds of all the other variables is greater, then it can + // never fall into a potential hall interval formed by the variable in the + // current window, so we can split the problem into independent parts. + if (lb <= min_lb + IntegerValue(num_in_window - 1)) { + ++num_in_window; + continue; + } + + // Process the current window. + if (num_in_window > 1) { + absl::Span window(&vars_[start], num_in_window); + if (!PropagateLowerBoundsInternal(min_lb, window)) { return false; } } - // Updates critical_intervals_. Note that we use the old lb, but that - // doesn't change the value of newly_covered. This block is what takes the - // most time. - int64 newly_covered; - const auto it = - critical_intervals_.GrowRightByOne(lb.value(), &newly_covered); - to_insert_.push_back({newly_covered, var}); - const IntegerValue end(it->end); + // Start of the next window. + start = i; + num_in_window = 1; + min_lb = lb; + } + + // Take care of the last window. + if (num_in_window > 1) { + absl::Span window(&vars_[start], num_in_window); + return PropagateLowerBoundsInternal(min_lb, window); + } + + return true; +} + +bool AllDifferentBoundsPropagator::PropagateLowerBoundsInternal( + IntegerValue min_lb, absl::Span vars) { + hall_starts_.clear(); + hall_ends_.clear(); + + // All cached lb in vars will be in [min_lb, min_lb + vars_.size()). + // Make sure we change our base_ so that GetIndex() fit in our buffers. + base_ = min_lb - IntegerValue(1); + + // Sparse cleaning of value_to_nodes_. + for (const int i : indices_to_clear_) { + index_to_var_[i] = kNoIntegerVariable; + } + indices_to_clear_.clear(); + + // Sort vars by increasing ub. + std::sort(vars.begin(), vars.end(), + [](VarValue a, VarValue b) { return a.ub < b.ub; }); + for (const VarValue entry : vars) { + const IntegerVariable var = entry.var; + + // Note that it is important to use the cache to make sure GetIndex() is + // not out of bound in case integer_trail_->LowerBound() changed when we + // pushed something. + const IntegerValue lb = entry.lb; + const int lb_index = GetIndex(lb); + const bool value_is_covered = PointIsPresent(lb_index); + + // Check if lb is in an Hall interval, and push it if this is the case. + if (value_is_covered) { + const int hall_index = + std::lower_bound(hall_ends_.begin(), hall_ends_.end(), lb) - + hall_ends_.begin(); + if (hall_index < hall_ends_.size() && hall_starts_[hall_index] <= lb) { + const IntegerValue hs = hall_starts_[hall_index]; + const IntegerValue he = hall_ends_[hall_index]; + FillHallReason(hs, he); + integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, hs)); + if (!integer_trail_->Enqueue( + IntegerLiteral::GreaterOrEqual(var, he + 1), + /*literal_reason=*/{}, integer_reason_)) { + return false; + } + } + } + + // Update our internal representation of the non-consecutive intervals. + // + // If lb is not used, we add a node there, otherwise we add it to the + // right of the interval that contains lb. In both cases, if there is an + // interval to the left (resp. right) we merge them. + int new_index = lb_index; + int start_index = lb_index; + int end_index = lb_index; + if (value_is_covered) { + start_index = FindStartIndexAndCompressPath(new_index); + new_index = index_to_end_index_[start_index] + 1; + end_index = new_index; + } else { + if (PointIsPresent(new_index - 1)) { + start_index = FindStartIndexAndCompressPath(new_index - 1); + } + } + if (PointIsPresent(new_index + 1)) { + end_index = index_to_end_index_[new_index + 1]; + index_to_start_index_[new_index + 1] = start_index; + } + + // Update the end of the representative. + index_to_end_index_[start_index] = end_index; + + // This is the only place where we "add" a new node. + { + index_to_start_index_[new_index] = start_index; + index_to_var_[new_index] = var; + indices_to_clear_.push_back(new_index); + } // We cannot have a conflict, because it should have beend detected before // by pushing an interval lower bound past its upper bound. + // + // TODO(user): Not 100% clear since pushing can have side-effect, maybe we + // should just report the conflict if it happens! + const IntegerValue end = GetValue(end_index); DCHECK_LE(end, integer_trail_->UpperBound(var)); // If we have a new Hall interval, add it to the set. Note that it will // always be last, and if it overlaps some previous Hall intervals, it // always overlaps them fully. - if (end == integer_trail_->UpperBound(var)) { - const IntegerValue start(it->start); + // + // Note: It is okay to not use entry.ub here if we want to fetch the last + // value, but in practice it shouldn't really change when we push a + // lower_bound and it is faster to use the cached entry. + if (end == entry.ub) { + const IntegerValue start = GetValue(start_index); while (!hall_starts_.empty() && start <= hall_starts_.back()) { hall_starts_.pop_back(); hall_ends_.pop_back(); @@ -511,8 +626,8 @@ bool AllDifferentBoundsPropagator::PropagateLowerBounds() { void AllDifferentBoundsPropagator::RegisterWith( GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); - for (const IntegerVariable& var : vars_) { - watcher->WatchIntegerVariable(var, id); + for (const VarValue entry : vars_) { + watcher->WatchIntegerVariable(entry.var, id); } watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); } diff --git a/ortools/sat/all_different.h b/ortools/sat/all_different.h index 25118abf52..66286addd4 100644 --- a/ortools/sat/all_different.h +++ b/ortools/sat/all_different.h @@ -24,7 +24,6 @@ #include "ortools/sat/integer.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" -#include "ortools/util/sorted_interval_list.h" namespace operations_research { namespace sat { @@ -140,14 +139,11 @@ class AllDifferentConstraint : PropagatorInterface { // deduce that the domain of the other variables cannot contains such Hall // interval. // -// We use a "simple" O(n log n) algorithm. +// We use a "fast" O(n log n) algorithm. // -// TODO(user): implement the faster algorithm described in: +// TODO(user): It might be difficult to find something faster than what is +// implemented here. Some related reference: // https://cs.uwaterloo.ca/~vanbeek/Publications/ijcai03_TR.pdf -// Note that the algorithms are similar, the gain comes by replacing our -// SortedDisjointIntervalList with a more customized class for our operations. -// It is even possible to get an O(n) complexity if the values of the bounds are -// in a range of size O(n). class AllDifferentBoundsPropagator : public PropagatorInterface { public: AllDifferentBoundsPropagator(const std::vector& vars, @@ -157,36 +153,71 @@ class AllDifferentBoundsPropagator : public PropagatorInterface { void RegisterWith(GenericLiteralWatcher* watcher); private: + // We locally cache the lb/ub for faster sorting and to guarantee some + // invariant when we push bounds. + struct VarValue { + IntegerVariable var; + IntegerValue lb; + IntegerValue ub; + }; + // Fills integer_reason_ with the reason why we have the given hall interval. void FillHallReason(IntegerValue hall_lb, IntegerValue hall_ub); - // Do half the job of Propagate(). + // Do half the job of Propagate(). This will split the variable into + // independent subset, and call PropagateLowerBoundsInternal() on each of + // them. bool PropagateLowerBounds(); + bool PropagateLowerBoundsInternal(IntegerValue min_lb, + absl::Span vars); + + // Internally, we will maintain a set of non-consecutive integer intervals of + // the form [start, end]. Each point (i.e. IntegerValue) of such interval will + // be associated to an unique variable and via an union-find algorithm point + // to its start. The end only make sense for representative. + // + // TODO(user): Because we don't use rank, we have a worst case complexity of + // O(n log n). We could try a normal Union-find data structure, but then we + // also have to maintain a start vector. + // + // Note that during the execution of the algorithm we start from empty + // intervals and finish with a set of points of size num_vars. + // + // The list of all points are maintained in the dense vectors index_to_*_ + // where we have remapped values to indices (with GetIndex()) to make sure it + // always fall into the correct range. + int FindStartIndexAndCompressPath(int index); + + int GetIndex(IntegerValue value) const { + DCHECK_GE(value, base_); + DCHECK_LT(value - base_, index_to_start_index_.size()); + return (value - base_).value(); + } + + IntegerValue GetValue(int index) const { return base_ + IntegerValue(index); } + + bool PointIsPresent(int index) const { + return index_to_var_[index] != kNoIntegerVariable; + } - std::vector vars_; - std::vector negated_vars_; IntegerTrail* integer_trail_; - // The sets of "critical" intervals. This has the same meaning as in the - // disjunctive constraint. - SortedDisjointIntervalList critical_intervals_; + // These vector will be either sorted by lb or by ub. + std::vector vars_; + std::vector negated_vars_; // The list of Hall intervalls detected so far, sorted. std::vector hall_starts_; std::vector hall_ends_; - // Members needed for explaining the propagation. - // - // The IntegerVariable in an hall interval [lb, ub] are the variables with key - // in [lb, ub] in this map. Note(user): if the set of bounds is small, we - // could use a std::vector here. The O(ub - lb) to create the reason is fine - // since this is the size of the reason. - // - // Optimization: we only insert the entry in the map lazily when the reason - // is needed. - int64 num_calls_; - std::vector> to_insert_; - absl::flat_hash_map value_to_variable_; + // Non-consecutive intervals related data-structures. + IntegerValue base_; + std::vector indices_to_clear_; + std::vector index_to_start_index_; + std::vector index_to_end_index_; + std::vector index_to_var_; + + // Temporary integer reason. std::vector integer_reason_; DISALLOW_COPY_AND_ASSIGN(AllDifferentBoundsPropagator); diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 50ccce5d8c..337c235f5b 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -61,6 +61,18 @@ NeighborhoodGeneratorHelper::NeighborhoodGeneratorHelper( } } } + + type_to_constraints_.clear(); + const int num_constraints = model_proto_.constraints_size(); + for (int c = 0; c < num_constraints; ++c) { + const int type = model_proto_.constraints(c).constraint_case(); + CHECK_GE(type, 0) << "Negative constraint case ??"; + CHECK_LT(type, 10000) << "Large constraint case ??"; + if (type >= type_to_constraints_.size()) { + type_to_constraints_.resize(type + 1); + } + type_to_constraints_[type].push_back(c); + } } bool NeighborhoodGeneratorHelper::IsActive(int var) const { @@ -111,20 +123,26 @@ CpModelProto NeighborhoodGeneratorHelper::RelaxGivenVariables( return FixGivenVariables(initial_solution, fixed_variables); } -CpModelProto SimpleNeighborhoodGenerator::Generate( - const CpSolverResponse& initial_solution, int64 seed, - double difficulty) const { +namespace { + +void GetRandomSubset(int seed, double relative_size, std::vector* base) { random_engine_t random; random.seed(seed); // TODO(user): we could generate this more efficiently than using random // shuffle. - std::vector fixed_variables = helper_.ActiveVariables(); + std::shuffle(base->begin(), base->end(), random); + const int target_size = std::ceil(relative_size * base->size()); + base->resize(target_size); +} - std::shuffle(fixed_variables.begin(), fixed_variables.end(), random); - const int target_size = - std::ceil((1.0 - difficulty) * fixed_variables.size()); - fixed_variables.resize(target_size); +} // namespace + +CpModelProto SimpleNeighborhoodGenerator::Generate( + const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const { + std::vector fixed_variables = helper_.ActiveVariables(); + GetRandomSubset(seed, 1.0 - difficulty, &fixed_variables); return helper_.FixGivenVariables(initial_solution, fixed_variables); } @@ -242,5 +260,91 @@ CpModelProto ConstraintGraphNeighborhoodGenerator::Generate( return helper_.RelaxGivenVariables(initial_solution, relaxed_variables); } +CpModelProto SchedulingNeighborhoodGenerator::Generate( + const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const { + std::set intervals_to_relax; + { + const auto span = helper_.TypeToConstraints(ConstraintProto::kInterval); + std::vector v(span.begin(), span.end()); + GetRandomSubset(seed, difficulty, &v); + intervals_to_relax.insert(v.begin(), v.end()); + } + + CpModelProto copy = helper_.ModelProto(); + + // Fix the presence/absence of non-relaxed intervals. + for (const int i : helper_.TypeToConstraints(ConstraintProto::kInterval)) { + if (intervals_to_relax.count(i)) continue; + + const ConstraintProto& interval_ct = copy.constraints(i); + if (interval_ct.enforcement_literal().empty()) continue; + + CHECK_EQ(interval_ct.enforcement_literal().size(), 1); + const int enforcement_ref = interval_ct.enforcement_literal(0); + const int enforcement_var = PositiveRef(enforcement_ref); + const int value = initial_solution.solution(enforcement_var); + + // Fix the value. + copy.mutable_variables(enforcement_var)->clear_domain(); + copy.mutable_variables(enforcement_var)->add_domain(value); + copy.mutable_variables(enforcement_var)->add_domain(value); + + // If the interval is ignored, skip for the loop below as there is no + // point adding precedence on it. + if (RefIsPositive(enforcement_ref) == (value == 0)) { + intervals_to_relax.insert(i); + } + } + + for (const int c : helper_.TypeToConstraints(ConstraintProto::kNoOverlap)) { + // Sort all non-relaxed intervals of this constraint by current start time. + std::vector> start_interval_pairs; + for (const int i : copy.constraints(c).no_overlap().intervals()) { + if (intervals_to_relax.count(i)) continue; + const ConstraintProto& interval_ct = copy.constraints(i); + + // TODO(user): we ignore size zero for now. + const int size_var = interval_ct.interval().size(); + if (initial_solution.solution(size_var) == 0) continue; + + const int start_var = interval_ct.interval().start(); + const int64 start_value = initial_solution.solution(start_var); + start_interval_pairs.push_back({start_value, i}); + } + std::sort(start_interval_pairs.begin(), start_interval_pairs.end()); + + // Add precedence between the remaining intervals, forcing their order. + for (int i = 0; i + 1 < start_interval_pairs.size(); ++i) { + const int before_var = + copy.constraints(start_interval_pairs[i].second).interval().end(); + const int after_var = copy.constraints(start_interval_pairs[i + 1].second) + .interval() + .start(); + CHECK_LE(initial_solution.solution(before_var), + initial_solution.solution(after_var)); + + LinearConstraintProto* linear = copy.add_constraints()->mutable_linear(); + linear->add_domain(kint64min); + linear->add_domain(0); + linear->add_vars(before_var); + linear->add_coeffs(1); + linear->add_vars(after_var); + linear->add_coeffs(-1); + } + } + + // Set the current solution as a hint. + // + // TODO(user): Move to common function? + copy.clear_solution_hint(); + for (int var = 0; var < copy.variables_size(); ++var) { + copy.mutable_solution_hint()->add_vars(var); + copy.mutable_solution_hint()->add_values(initial_solution.solution(var)); + } + + return copy; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index cfa91a2cb3..068a669abd 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -16,6 +16,7 @@ #include +#include "absl/types/span.h" #include "ortools/base/integral_types.h" #include "ortools/sat/cp_model.pb.h" @@ -60,6 +61,13 @@ class NeighborhoodGeneratorHelper { return var_to_constraint_; } + // Returns all the constraints indices of a given type. + const absl::Span TypeToConstraints( + ConstraintProto::ConstraintCase type) const { + if (type >= type_to_constraints_.size()) return {}; + return absl::MakeSpan(type_to_constraints_[type]); + } + // The initial problem. const CpModelProto& ModelProto() const { return model_proto_; } @@ -69,6 +77,9 @@ class NeighborhoodGeneratorHelper { const CpModelProto& model_proto_; + // Constraints by types. + std::vector> type_to_constraints_; + // Variable-Constraint graph. std::vector> constraint_to_var_; std::vector> var_to_constraint_; @@ -150,6 +161,21 @@ class ConstraintGraphNeighborhoodGenerator : public NeighborhoodGenerator { double difficulty) const final; }; +// Only make sense for scheduling problem. This select a random set of interval +// of the problem according to the difficulty. Then, for each no_overlap +// constraints, it adds strict relation order between the non-relaxed intervals. +// +// TODO(user): Also deal with cumulative constraint. +class SchedulingNeighborhoodGenerator : public NeighborhoodGenerator { + public: + explicit SchedulingNeighborhoodGenerator( + NeighborhoodGeneratorHelper const* helper, const std::string& name) + : NeighborhoodGenerator(name, helper) {} + + CpModelProto Generate(const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const final; +}; + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index bbe6bb9636..b8b2cab0f4 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -497,9 +497,9 @@ std::string CpSolverResponseStats(const CpSolverResponse& response) { absl::StrAppend(&result, "\nobjective: NA"); absl::StrAppend(&result, "\nbest_bound: NA"); } else { - absl::StrAppendFormat(&result, "\nobjective: %g", + absl::StrAppendFormat(&result, "\nobjective: %.9g", response.objective_value()); - absl::StrAppendFormat(&result, "\nbest_bound: %g", + absl::StrAppendFormat(&result, "\nbest_bound: %.9g", response.best_objective_bound()); } @@ -513,10 +513,10 @@ std::string CpSolverResponseStats(const CpSolverResponse& response) { "\npropagations: ", response.num_binary_propagations()); absl::StrAppend( &result, "\ninteger_propagations: ", response.num_integer_propagations()); - absl::StrAppendFormat(&result, "\nwalltime: %g", response.wall_time()); - absl::StrAppendFormat(&result, "\nusertime: %g", response.user_time()); - absl::StrAppendFormat(&result, "\ndeterministic_time: %g", - response.deterministic_time()); + absl::StrAppend(&result, "\nwalltime: ", response.wall_time()); + absl::StrAppend(&result, "\nusertime: ", response.user_time()); + absl::StrAppend(&result, + "\ndeterministic_time: ", response.deterministic_time()); absl::StrAppend(&result, "\n"); return result; } @@ -1900,8 +1900,16 @@ CpSolverResponse SolveCpModelWithLNS( const bool focus_on_decision_variables = parameters->lns_focus_on_decision_variables(); + // TODO(user): Find a way to propagate the level zero bounds from the other + // worker inside this base LNS problem. + const NeighborhoodGeneratorHelper helper(&model_proto, + focus_on_decision_variables); + // For now we will just alternate between our possible neighborhoods. - NeighborhoodGeneratorHelper helper(&model_proto, focus_on_decision_variables); + // + // TODO(user): select with higher probability the one that seems to work best? + // We can evaluate this compared to the best solution at the time of the + // neighborhood creation. std::vector> generators; generators.push_back( absl::make_unique(&helper, "rnd_lns")); @@ -1909,9 +1917,17 @@ CpSolverResponse SolveCpModelWithLNS( &helper, "var_lns")); generators.push_back(absl::make_unique( &helper, "cst_lns")); + if (!helper.TypeToConstraints(ConstraintProto::kNoOverlap).empty()) { + generators.push_back(absl::make_unique( + &helper, "scheduling_lns")); + } // The "optimal" difficulties do not have to be the same for different - // generators. TODO(user): move this inside the generator API? + // generators. + // + // TODO(user): This should be shared across LNS threads! create a thread + // safe class that accept signals of the form (difficulty and result) and + // properly update the difficulties. std::vector difficulties(generators.size(), AdaptiveParameterValue(0.5)); diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index 7f4b258581..62fee7ff69 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -38,18 +38,6 @@ NonOverlappingRectanglesPropagator::NonOverlappingRectanglesPropagator( strict_(strict), integer_trail_(integer_trail) { CHECK_GT(num_boxes_, 0); - neighbors_.resize(num_boxes_ * (num_boxes_ - 1)); - neighbors_begins_.resize(num_boxes_); - neighbors_ends_.resize(num_boxes_); - for (int i = 0; i < num_boxes_; ++i) { - const int begin = i * (num_boxes_ - 1); - neighbors_begins_[i] = begin; - neighbors_ends_[i] = begin + (num_boxes_ - 1); - for (int j = 0; j < num_boxes_; ++j) { - if (j == i) continue; - neighbors_[begin + (j > i ? j - 1 : j)] = j; - } - } } NonOverlappingRectanglesPropagator::~NonOverlappingRectanglesPropagator() {} @@ -68,7 +56,7 @@ bool NonOverlappingRectanglesPropagator::Propagate() { y_.SetTimeDirection(true); for (int box = 0; box < num_boxes_; ++box) { - if (!strict_ && cached_areas_[box] == 0) continue; + if (cached_areas_[box] == 0) continue; RETURN_IF_FALSE(FailWhenEnergyIsTooLarge(box)); } @@ -80,9 +68,6 @@ void NonOverlappingRectanglesPropagator::RegisterWith( const int id = watcher->Register(this); x_.WatchAllTasks(id, watcher); y_.WatchAllTasks(id, watcher); - for (int t = 0; t < num_boxes_; ++t) { - watcher->RegisterReversibleInt(id, &neighbors_ends_[t]); - } watcher->SetPropagatorPriority(id, 3); } @@ -106,57 +91,40 @@ IntegerValue DistanceToBoundingInterval(IntegerValue min_a, IntegerValue max_a, } // namespace -void NonOverlappingRectanglesPropagator::UpdateNeighbors(int box) { - tmp_removed_.clear(); +void NonOverlappingRectanglesPropagator::SortNeighbors(int box) { cached_distance_to_bounding_box_.resize(num_boxes_); + neighbors_.clear(); const IntegerValue box_x_min = x_.StartMin(box); const IntegerValue box_x_max = x_.EndMax(box); const IntegerValue box_y_min = y_.StartMin(box); const IntegerValue box_y_max = y_.EndMax(box); - int new_index = neighbors_begins_[box]; - const int end = neighbors_ends_[box]; - for (int i = new_index; i < end; ++i) { - const int other = neighbors_[i]; + + for (int other = 0; other < num_boxes_; ++other) { + if (other == box) continue; + if (cached_areas_[other] == 0) continue; const IntegerValue other_x_min = x_.StartMin(other); const IntegerValue other_x_max = x_.EndMax(other); - if (IntervalAreDisjointForSure(box_x_min, box_x_max, other_x_min, - other_x_max)) { - tmp_removed_.push_back(other); - continue; - } - const IntegerValue other_y_min = y_.StartMin(other); const IntegerValue other_y_max = y_.EndMax(other); - if (IntervalAreDisjointForSure(box_y_min, box_y_max, other_y_min, - other_y_max)) { - tmp_removed_.push_back(other); - continue; - } - neighbors_[new_index++] = other; + neighbors_.push_back(other); cached_distance_to_bounding_box_[other] = std::max(DistanceToBoundingInterval(box_x_min, box_x_max, other_x_min, other_x_max), DistanceToBoundingInterval(box_y_min, box_y_max, other_y_min, other_y_max)); } - neighbors_ends_[box] = new_index; - for (int i = 0; i < tmp_removed_.size();) { - neighbors_[new_index++] = tmp_removed_[i++]; - } - IncrementalSort(neighbors_.begin() + neighbors_begins_[box], - neighbors_.begin() + neighbors_ends_[box], - [this](int i, int j) { - return cached_distance_to_bounding_box_[i] < - cached_distance_to_bounding_box_[j]; - }); + IncrementalSort(neighbors_.begin(), neighbors_.begin(), [this](int i, int j) { + return cached_distance_to_bounding_box_[i] < + cached_distance_to_bounding_box_[j]; + }); } bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) { // Note that we only consider the smallest dimension of each boxes here. - UpdateNeighbors(box); + SortNeighbors(box); IntegerValue area_min_x = x_.StartMin(box); IntegerValue area_max_x = x_.EndMax(box); IntegerValue area_min_y = y_.StartMin(box); @@ -165,9 +133,7 @@ bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) { IntegerValue sum_of_areas = cached_areas_[box]; IntegerValue total_sum_of_areas = sum_of_areas; - const int end = neighbors_ends_[box]; - for (int i = neighbors_begins_[box]; i < end; ++i) { - const int other = neighbors_[i]; + for (const int other : neighbors_) { total_sum_of_areas += cached_areas_[other]; } @@ -182,9 +148,9 @@ bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) { // TODO(user): Is there a better order, maybe sort by distance // with the current box. - for (int i = neighbors_begins_[box]; i < end; ++i) { + for (int i = 0; i < neighbors_.size(); ++i) { const int other = neighbors_[i]; - if (cached_areas_[other] == 0) continue; + CHECK_GT(cached_areas_[other], 0); // Update Bounding box. area_min_x = std::min(area_min_x, x_.StartMin(other)); @@ -205,10 +171,8 @@ bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) { x_.ClearReason(); y_.ClearReason(); add_box_energy_in_rectangle_reason(box); - for (int j = neighbors_begins_[box]; j <= i; ++j) { - const int other = neighbors_[j]; - if (cached_areas_[other] == 0) continue; - add_box_energy_in_rectangle_reason(other); + for (int j = 0; j <= i; ++j) { + add_box_energy_in_rectangle_reason(neighbors_[j]); } x_.ImportOtherReasons(y_); return x_.ReportConflict(); @@ -290,7 +254,6 @@ bool CheckOverload(bool time_direction, IntegerValue other_time, other->AddEndMinReason(task, other_time + 1); } } - // LOG(INFO) << "Overload: " << num_events - critical_event; helper->ImportOtherReasons(*other); return helper->ReportConflict(); } @@ -298,12 +261,30 @@ bool CheckOverload(bool time_direction, IntegerValue other_time, return true; } +void AddOtherReasons(const std::vector& tasks, IntegerValue other_time, + int main_task, SchedulingConstraintHelper* other) { + other->ClearReason(); + bool main_task_seen = false; + for (const int task : tasks) { + other->AddStartMaxReason(task, other_time); + other->AddEndMinReason(task, other_time + 1); + if (task == main_task) { + main_task_seen = true; + } + } + if (!main_task_seen) { + other->AddStartMaxReason(main_task, other_time); + other->AddEndMinReason(main_task, other_time + 1); + } +} + bool DetectPrecedences(bool time_direction, IntegerValue other_time, - const absl::flat_hash_set& boxes, + const absl::flat_hash_set& active_boxes, SchedulingConstraintHelper* helper, SchedulingConstraintHelper* other) { helper->SetTimeDirection(time_direction); other->SetTimeDirection(true); + const auto& task_by_increasing_end_min = helper->TaskByIncreasingEndMin(); const auto& task_by_decreasing_start_max = helper->TaskByDecreasingStartMax(); @@ -314,7 +295,7 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time, const int t = task_time.task_index; const IntegerValue end_min = task_time.time; - if (helper->IsAbsent(t) || !gtl::ContainsKey(boxes, t)) continue; + if (helper->IsAbsent(t) || !gtl::ContainsKey(active_boxes, t)) continue; while (queue_index >= 0) { const auto to_insert = task_by_decreasing_start_max[queue_index]; @@ -322,7 +303,7 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time, const IntegerValue start_max = to_insert.time; if (end_min <= start_max) break; - if (gtl::ContainsKey(boxes, task_index)) { + if (gtl::ContainsKey(active_boxes, task_index)) { CHECK(helper->IsPresent(task_index)); task_set.AddEntry({task_index, helper->ShiftedStartMin(task_index), helper->DurationMin(task_index)}); @@ -343,44 +324,22 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time, if (end_min_of_critical_tasks > helper->StartMin(t)) { const std::vector& sorted_tasks = task_set.SortedTasks(); helper->ClearReason(); - other->ClearReason(); // We need: // - StartMax(ct) < EndMin(t) for the detectable precedence. // - StartMin(ct) > window_start for the end_min_of_critical_tasks reason. const IntegerValue window_start = sorted_tasks[critical_index].start_min; - const int num_conflicting_tasks = sorted_tasks.size() - critical_index; - if (num_conflicting_tasks < 4) { - for (int i = critical_index; i < sorted_tasks.size(); ++i) { - const int ct = sorted_tasks[i].task; - if (ct == t) continue; - other->AddReasonForBeingBefore(t, ct); - other->AddReasonForBeingBefore(ct, t); - } - for (int i = critical_index; i + 1 < sorted_tasks.size(); ++i) { - const int ct1 = sorted_tasks[i].task; - if (ct1 == t) continue; - for (int j = i + 1; j < sorted_tasks.size(); ++j) { - const int ct2 = sorted_tasks[j].task; - if (ct2 == t) continue; - other->AddReasonForBeingBefore(ct1, ct2); - other->AddReasonForBeingBefore(ct2, ct1); - } - } - } else { - for (int i = critical_index; i < sorted_tasks.size(); ++i) { - const int ct = sorted_tasks[i].task; - if (ct == t) continue; - other->AddStartMaxReason(ct, other_time); - other->AddEndMinReason(ct, other_time + 1); - } - other->AddStartMaxReason(t, other_time); - other->AddEndMinReason(t, other_time + 1); + std::vector tasks(sorted_tasks.size() - critical_index); + for (int i = critical_index; i < sorted_tasks.size(); ++i) { + tasks[i - critical_index] = sorted_tasks[i].task; } + + AddOtherReasons(tasks, other_time, t, other); + for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; if (ct == t) continue; - CHECK(gtl::ContainsKey(boxes, ct)); + CHECK(gtl::ContainsKey(active_boxes, ct)); CHECK(helper->IsPresent(ct)); helper->AddEnergyAfterReason(ct, sorted_tasks[i].duration_min, window_start); @@ -390,7 +349,7 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time, // Add the reason for t (we only need the end-min). helper->AddEndMinReason(t, end_min); - // Import other reasons. + // Import reasons from the 'other' dimension. helper->ImportOtherReasons(*other); // This augment the start-min of t and subsequently it can augment the @@ -409,7 +368,6 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time, } return true; } - } // namespace bool NonOverlappingRectanglesPropagator::PropagateOnProjections() { @@ -472,7 +430,26 @@ bool NonOverlappingRectanglesPropagator::FindMandatoryBoxesOnOneDimension( } for (const auto& it : mandatory_boxes) { - RETURN_IF_FALSE(PropagateMandatoryBoxesOnOneDimension(it.first, it.second, + // Compute the 'canonical' line to use when explaining that boxes overlap + // on the 'other' dimension. + IntegerValue lb(kint64min); + IntegerValue ub(kint64max); + for (const int task : it.second) { + lb = std::max(lb, other->StartMax(task)); + ub = std::min(ub, other->EndMin(task)); + } + + const IntegerValue span = ub - lb + 1; + IntegerValue selected = lb; + for (int shift = 30; shift >= 0; --shift) { + const IntegerValue mask(static_cast(1) >> shift); + if (mask <= span) { + selected = (ub / mask) * mask; + break; + } + } + // And propagate. + RETURN_IF_FALSE(PropagateMandatoryBoxesOnOneDimension(selected, it.second, helper, other)); } return true; @@ -485,109 +462,6 @@ bool NonOverlappingRectanglesPropagator::PropagateMandatoryBoxesOnOneDimension( RETURN_IF_FALSE(CheckOverload(true, event, active_boxes, helper, other)); RETURN_IF_FALSE(DetectPrecedences(true, event, active_boxes, helper, other)); RETURN_IF_FALSE(DetectPrecedences(false, event, active_boxes, helper, other)); - // RETURN_IF_FALSE(NaivePush(true, event, active_boxes, helper, other)); - // RETURN_IF_FALSE(NaivePush(false, event, active_boxes, helper, other)); - return true; -} - -bool NonOverlappingRectanglesPropagator::NaivePush( - bool time_direction, IntegerValue other_time, - const absl::flat_hash_set& boxes, SchedulingConstraintHelper* helper, - SchedulingConstraintHelper* other) { - // const auto left_box_before_right_box = - // [](int left, int right, SchedulingConstraintHelper* helper) { - // // left box pushes right box. - // const IntegerValue left_end_min = helper->EndMin(left); - // if (left_end_min > helper->StartMin(right)) { - // // Store reasons state. - // const int literal_size = helper->MutableLiteralReason()->size(); - // const int integer_size = helper->MutableIntegerReason()->size(); - // - // helper->AddEndMinReason(left, left_end_min); - // if (!helper->IncreaseStartMin(right, left_end_min)) { - // return false; - // } - // - // // Restore the reasons to the state before the increase of the - // start. helper->MutableLiteralReason()->resize(literal_size); - // helper->MutableIntegerReason()->resize(integer_size); - // } - // - // // right box pushes left box. - // const IntegerValue right_start_max = helper->StartMax(right); - // if (right_start_max < helper->EndMax(left)) { - // helper->AddStartMaxReason(right, right_start_max); - // return helper->DecreaseEndMax(left, right_start_max); - // } - // - // return true; - // }; - // - // helper->SetTimeDirection(time_direction); - // other->SetTimeDirection(true); - // const auto& task_by_increasing_end_min = - // helper->TaskByIncreasingEndMin(); const auto& - // task_by_decreasing_start_max = - // helper->TaskByDecreasingStartMax(); - // - // int queue_index = helper->NumTasks() - 1; - // for (const auto task_time : task_by_increasing_end_min) { - // const int t = task_time.task_index; - // const IntegerValue end_min = task_time.time; - // - // if (!gtl::ContainsKey(boxes, t)) continue; - // - // - // for (int b1 = 0; b1 + 1 < boxes.size(); ++b1) { - // const int box1 = boxes[b1]; - // if (!strict_ && cached_areas_[box1] == 0) continue; - // for (int b2 = b1 + 1; b2 < boxes.size(); ++b2) { - // const int box2 = boxes[b2]; - // if (!strict_ && cached_areas_[box2] == 0) continue; - // // For each direction and each order, we test if the boxes can be - // // disjoint. - // const int state = (helper->EndMin(box1) <= helper->StartMax(box2)) + - // 2 * (helper->EndMin(box2) <= - // helper->StartMax(box1)); - // - // if (state == 3) continue; - // // Clean up reasons. - // helper->ClearReason(); - // other->ClearReason(); - // - // // This is an "hack" to be able to easily test for none or for one - // // and only one of the conditions below. - // switch (state) { - // case 0: { - // helper->AddReasonForBeingBefore(box1, box2); - // helper->AddReasonForBeingBefore(box2, box1); - // other->AddReasonForBeingBefore(box1, box2); - // other->AddReasonForBeingBefore(box2, box1); - // helper->ImportOtherReasons(*other); - // return helper->ReportConflict(); - // } - // case 1: { - // other->AddReasonForBeingBefore(box1, box2); - // other->AddReasonForBeingBefore(box2, box1); - // helper->AddReasonForBeingBefore(box1, box2); - // helper->ImportOtherReasons(*other); - // RETURN_IF_FALSE(left_box_before_right_box(box1, box2, helper)); - // break; - // } - // case 2: { - // other->AddReasonForBeingBefore(box1, box2); - // other->AddReasonForBeingBefore(box2, box1); - // helper->AddReasonForBeingBefore(box2, box1); - // helper->ImportOtherReasons(*other); - // RETURN_IF_FALSE(left_box_before_right_box(box2, box1, helper)); - // break; - // } - // default: { - // break; - // } - // } - // } - // } return true; } diff --git a/ortools/sat/diffn.h b/ortools/sat/diffn.h index 6bdfd20b25..a00131c845 100644 --- a/ortools/sat/diffn.h +++ b/ortools/sat/diffn.h @@ -52,11 +52,7 @@ class NonOverlappingRectanglesPropagator : public PropagatorInterface { const std::vector& boxes, SchedulingConstraintHelper* active, SchedulingConstraintHelper* other); - bool NaivePush(bool time_direction, IntegerValue other_time, - const absl::flat_hash_set& boxes, - SchedulingConstraintHelper* helper, - SchedulingConstraintHelper* other_helper); - void UpdateNeighbors(int box); + void SortNeighbors(int box); bool FailWhenEnergyIsTooLarge(int box); const int num_boxes_; @@ -65,13 +61,7 @@ class NonOverlappingRectanglesPropagator : public PropagatorInterface { const bool strict_; IntegerTrail* integer_trail_; - // The neighbors_ of a box will be in - // [neighbors_[begin[box]], neighbors_[end[box]]) std::vector neighbors_; - std::vector neighbors_begins_; - std::vector neighbors_ends_; - std::vector tmp_removed_; - std::vector cached_areas_; std::vector cached_distance_to_bounding_box_; diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index 164bbaec9c..aac7395c65 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -25,25 +25,94 @@ namespace operations_research { namespace sat { +LiteralIndex AtMinValue(IntegerVariable var, Model* model) { + IntegerEncoder* const integer_encoder = model->GetOrCreate(); + IntegerTrail* const integer_trail = model->GetOrCreate(); + + DCHECK(!integer_trail->IsCurrentlyIgnored(var)); + const IntegerValue lb = integer_trail->LowerBound(var); + DCHECK_LE(lb, integer_trail->UpperBound(var)); + if (lb == integer_trail->UpperBound(var)) return kNoLiteralIndex; + return integer_encoder + ->GetOrCreateAssociatedLiteral(IntegerLiteral::LowerOrEqual(var, lb)) + .Index(); +} + +LiteralIndex GreaterOrEqualToMiddleValue(IntegerVariable var, Model* model) { + auto* encoder = model->GetOrCreate(); + auto* trail = model->GetOrCreate(); + auto* integer_trail = model->GetOrCreate(); + const IntegerValue var_lb = integer_trail->LowerBound(var); + const IntegerValue var_ub = integer_trail->UpperBound(var); + CHECK_LT(var_lb, var_ub); + + const IntegerValue chosen_value = + var_lb + std::max(IntegerValue(1), (var_ub - var_lb) / IntegerValue(2)); + const Literal ge = encoder->GetOrCreateAssociatedLiteral( + IntegerLiteral::GreaterOrEqual(var, chosen_value)); + CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable())); + VLOG(2) << "Chosen " << var << " >= " << chosen_value; + return ge.Index(); +} + +LiteralIndex SplitDomainUsingLpValue(IntegerVariable var, Model* model) { + auto* encoder = model->GetOrCreate(); + auto* trail = model->GetOrCreate(); + auto* integer_trail = model->GetOrCreate(); + auto* lp_dispatcher = model->GetOrCreate(); + + const IntegerVariable positive_var = + VariableIsPositive(var) ? var : NegationOf(var); + DCHECK(!integer_trail->IsCurrentlyIgnored(positive_var)); + LinearProgrammingConstraint* lp = + gtl::FindWithDefault(*lp_dispatcher, positive_var, nullptr); + if (lp == nullptr) return kNoLiteralIndex; + const IntegerValue value = IntegerValue( + static_cast(std::round(lp->GetSolutionValue(positive_var)))); + + // Because our lp solution might be from higher up in the tree, it + // is possible that value is now outside the domain of positive_var. + // In this case, we just revert to the current literal. + const IntegerValue lb = integer_trail->LowerBound(positive_var); + const IntegerValue ub = integer_trail->UpperBound(positive_var); + + // We try first (<= value), but if this do not reduce the domain we + // try to enqueue (>= value). Note that even for domain with hole, + // since we know that this variable is not fixed, one of the two + // alternative must reduce the domain. + // + // TODO(user): use GetOrCreateLiteralAssociatedToEquality() instead? + // It may replace two decision by only one. However this function + // cannot currently be called during search, but that should be easy + // enough to fix. + if (value >= lb && value < ub) { + const Literal le = encoder->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(positive_var, value)); + CHECK(!trail->Assignment().VariableIsAssigned(le.Variable())); + return le.Index(); + } + if (value > lb && value <= ub) { + const Literal ge = encoder->GetOrCreateAssociatedLiteral( + IntegerLiteral::GreaterOrEqual(positive_var, value)); + CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable())); + return ge.Index(); + } + return kNoLiteralIndex; +} + // TODO(user): the complexity caused by the linear scan in this heuristic and // the one below is ok when search_branching is set to SAT_SEARCH because it is // not executed often, but otherwise it is done for each search decision, // which seems expensive. Improve. std::function FirstUnassignedVarAtItsMinHeuristic( const std::vector& vars, Model* model) { - IntegerEncoder* const integer_encoder = model->GetOrCreate(); - IntegerTrail* const integer_trail = model->GetOrCreate(); - return [integer_encoder, integer_trail, /*copy*/ vars]() { + return [/*copy*/ vars, model]() { + IntegerTrail* const integer_trail = model->GetOrCreate(); for (const IntegerVariable var : vars) { // Note that there is no point trying to fix a currently ignored variable. if (integer_trail->IsCurrentlyIgnored(var)) continue; - const IntegerValue lb = integer_trail->LowerBound(var); - if (lb < integer_trail->UpperBound(var)) { - return integer_encoder - ->GetOrCreateAssociatedLiteral( - IntegerLiteral::LowerOrEqual(var, lb)) - .Index(); - } + const LiteralIndex decision = AtMinValue(var, model); + if (decision != kNoLiteralIndex) return decision; } return kNoLiteralIndex; }; @@ -51,9 +120,8 @@ std::function FirstUnassignedVarAtItsMinHeuristic( std::function UnassignedVarWithLowestMinAtItsMinHeuristic( const std::vector& vars, Model* model) { - IntegerEncoder* const integer_encoder = model->GetOrCreate(); - IntegerTrail* const integer_trail = model->GetOrCreate(); - return [integer_encoder, integer_trail, /*copy */ vars]() { + return [/*copy */ vars, model]() { + IntegerTrail* const integer_trail = model->GetOrCreate(); IntegerVariable candidate = kNoIntegerVariable; IntegerValue candidate_lb; for (const IntegerVariable var : vars) { @@ -66,10 +134,7 @@ std::function UnassignedVarWithLowestMinAtItsMinHeuristic( } } if (candidate == kNoIntegerVariable) return kNoLiteralIndex; - return integer_encoder - ->GetOrCreateAssociatedLiteral( - IntegerLiteral::LowerOrEqual(candidate, candidate_lb)) - .Index(); + return AtMinValue(candidate, model); }; } @@ -96,10 +161,6 @@ std::function SatSolverHeuristic(Model* model) { } std::function PseudoCost(Model* model) { - auto* encoder = model->GetOrCreate(); - auto* trail = model->GetOrCreate(); - auto* integer_trail = model->GetOrCreate(); - const ObjectiveSynchronizationHelper* helper = model->Get(); const bool has_objective = @@ -111,29 +172,14 @@ std::function PseudoCost(Model* model) { PseudoCosts* pseudo_costs = model->GetOrCreate(); - // NOTE: The algorithm to choose the value for a variable is kept outside of - // PseudoCosts class since this is an independent heuristic. For now this is - // only used for pseudo costs however this can be refactored a bit for other - // branching strategies to use. return [=]() { const IntegerVariable chosen_var = pseudo_costs->GetBestDecisionVar(); if (chosen_var == kNoIntegerVariable) return kNoLiteralIndex; - const IntegerValue chosen_var_lb = integer_trail->LowerBound(chosen_var); - const IntegerValue chosen_var_ub = integer_trail->UpperBound(chosen_var); - CHECK_LT(chosen_var_lb, chosen_var_ub); // TODO(user): Experiment with different heuristics for choosing // value. - const IntegerValue chosen_value = - chosen_var_lb + - std::max(IntegerValue(1), - (chosen_var_ub - chosen_var_lb) / IntegerValue(2)); - const Literal ge = encoder->GetOrCreateAssociatedLiteral( - IntegerLiteral::GreaterOrEqual(chosen_var, chosen_value)); - CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable())); - VLOG(2) << "Chosen " << chosen_var << " >= " << chosen_value; - return ge.Index(); + return GreaterOrEqualToMiddleValue(chosen_var, model); }; } @@ -214,16 +260,50 @@ std::function FollowHint( }; } -std::function ExploitLpSolution( - std::function heuristic, Model* model) { +LiteralIndex ExploitLpSolution(const LiteralIndex decision, Model* model) { auto* encoder = model->GetOrCreate(); - auto* trail = model->GetOrCreate(); auto* integer_trail = model->GetOrCreate(); - auto* lp_dispatcher = model->GetOrCreate(); + auto* lp_constraints = model->GetOrCreate(); const SatParameters& parameters = *(model->GetOrCreate()); + if (decision == kNoLiteralIndex) { + return decision; + } + + bool lp_solution_is_exploitable = true; + // TODO(user,user): When we have more than one LP, their set of variable + // is always disjoint. So we could still change the polarity if the next + // variable we branch on is part of a LP that has a solution. + for (LinearProgrammingConstraint* lp : *lp_constraints) { + if (!lp->HasSolution() || + !(parameters.exploit_all_lp_solution() || lp->SolutionIsInteger())) { + lp_solution_is_exploitable = false; + break; + } + } + if (lp_solution_is_exploitable) { + for (const IntegerLiteral l : + encoder->GetAllIntegerLiterals(Literal(decision))) { + const IntegerVariable positive_var = + VariableIsPositive(l.var) ? l.var : NegationOf(l.var); + if (integer_trail->IsCurrentlyIgnored(positive_var)) continue; + const LiteralIndex lp_decision = + SplitDomainUsingLpValue(positive_var, model); + if (lp_decision != kNoLiteralIndex) { + return lp_decision; + } + } + } + return decision; +} + +std::function ExploitLpSolution( + std::function heuristic, Model* model) { + auto* lp_constraints = + model->GetOrCreate(); + // Use the normal heuristic if the LP(s) do not seem to cover enough variables // to be relevant. // TODO(user): Instead, try and depending on the result call it again or not? @@ -237,86 +317,9 @@ std::function ExploitLpSolution( if (num_integer_variables > 2 * num_lp_variables) return heuristic; } - bool last_decision_followed_lp = false; - int old_level = 0; - double old_obj = 0.0; return [=]() mutable { const LiteralIndex decision = heuristic(); - if (decision == kNoLiteralIndex) { - if (last_decision_followed_lp) { - VLOG(1) << "Integer LP solution is feasible, level:" << old_level - << "->" << trail->CurrentDecisionLevel() << " obj:" << old_obj; - } - last_decision_followed_lp = false; - return kNoLiteralIndex; - } - - bool lp_solution_is_exploitable = true; - double obj = 0.0; - // TODO(user,user): When we have more than one LP, their set of variable - // is always disjoint. So we could still change the polarity if the next - // variable we branch on is part of a LP that has a solution. - for (LinearProgrammingConstraint* lp : *lp_constraints) { - if (!lp->HasSolution() || - !(parameters.exploit_all_lp_solution() || lp->SolutionIsInteger())) { - lp_solution_is_exploitable = false; - break; - } - obj += lp->SolutionObjectiveValue(); - } - if (lp_solution_is_exploitable) { - if (!last_decision_followed_lp || obj != old_obj) { - old_level = trail->CurrentDecisionLevel(); - old_obj = obj; - VLOG(2) << "Integer LP solution at level:" << old_level - << " obj:" << old_obj; - } - for (const IntegerLiteral l : - encoder->GetAllIntegerLiterals(Literal(decision))) { - const IntegerVariable positive_var = - VariableIsPositive(l.var) ? l.var : NegationOf(l.var); - if (integer_trail->IsCurrentlyIgnored(positive_var)) continue; - LinearProgrammingConstraint* lp = - gtl::FindWithDefault(*lp_dispatcher, positive_var, nullptr); - if (lp != nullptr) { - const IntegerValue value = IntegerValue(static_cast( - std::round(lp->GetSolutionValue(positive_var)))); - - // Because our lp solution might be from higher up in the tree, it - // is possible that value is now outside the domain of positive_var. - // In this case, we just revert to the current literal. - const IntegerValue lb = integer_trail->LowerBound(positive_var); - const IntegerValue ub = integer_trail->UpperBound(positive_var); - - // We try first (<= value), but if this do not reduce the domain we - // try to enqueue (>= value). Note that even for domain with hole, - // since we know that this variable is not fixed, one of the two - // alternative must reduce the domain. - // - // TODO(user): use GetOrCreateLiteralAssociatedToEquality() instead? - // It may replace two decision by only one. However this function - // cannot currently be called during search, but that should be easy - // enough to fix. - if (value >= lb && value < ub) { - const Literal le = encoder->GetOrCreateAssociatedLiteral( - IntegerLiteral::LowerOrEqual(positive_var, value)); - CHECK(!trail->Assignment().VariableIsAssigned(le.Variable())); - last_decision_followed_lp = true; - return le.Index(); - } - if (value > lb && value <= ub) { - const Literal ge = encoder->GetOrCreateAssociatedLiteral( - IntegerLiteral::GreaterOrEqual(positive_var, value)); - CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable())); - last_decision_followed_lp = true; - return ge.Index(); - } - } - } - } - - last_decision_followed_lp = false; - return decision; + return ExploitLpSolution(decision, model); }; } @@ -598,7 +601,7 @@ SatSolver::Status SolveIntegerProblemWithLazyEncoding(Model* model) { void LogNewSolution(const std::string& event_or_solution_count, double time_in_seconds, double obj_lb, double obj_ub, const std::string& solution_info) { - LOG(INFO) << absl::StrFormat("#%-5s %6.2fs obj:[%g,%g] %s", + LOG(INFO) << absl::StrFormat("#%-5s %6.2fs obj:[%.9g,%.9g] %s", event_or_solution_count, time_in_seconds, obj_lb, obj_ub, solution_info); } diff --git a/ortools/sat/integer_search.h b/ortools/sat/integer_search.h index e6b36986fd..ce37d24f74 100644 --- a/ortools/sat/integer_search.h +++ b/ortools/sat/integer_search.h @@ -17,11 +17,24 @@ #include #include "ortools/sat/integer.h" +#include "ortools/sat/sat_base.h" #include "ortools/sat/sat_solver.h" namespace operations_research { namespace sat { +// Returns decision corresponding to var at its lower bound. Returns +// kNoLiteralIndex if the variable is fixed. +LiteralIndex AtMinValue(IntegerVariable var, Model* model); + +// Returns decision corresponding to var <= round(lp_value). If the variable +// does not appear in the LP, this method returns kNoLiteralIndex. +LiteralIndex SplitDomainUsingLpValue(IntegerVariable var, Model* model); + +// Returns decision corresponding to var >= lb + max(1, (ub - lb) / 2). It also +// CHECKs that the variable is not fixed. +LiteralIndex GreaterOrEqualToMiddleValue(IntegerVariable var, Model* model); + // Decision heuristic for SolveIntegerProblemWithLazyEncoding(). Returns a // function that will return the literal corresponding to the fact that the // first currently non-fixed variable value is <= its min. The function will @@ -75,6 +88,11 @@ std::function PseudoCost(Model* model); std::function ExploitLpSolution( std::function heuristic, Model* model); +// Similar to ExploitLpSolution(). Takes LiteralIndex as base decision and +// changes change the returned decision to AtLpValue() of the underlying integer +// variable if LP solution is exploitable. +LiteralIndex ExploitLpSolution(const LiteralIndex decision, Model* model); + // A restart policy that restarts every k failures. std::function RestartEveryKFailures(int k, SatSolver* solver); diff --git a/ortools/util/sorted_interval_list.cc b/ortools/util/sorted_interval_list.cc index eeb50980e9..da67c6cf00 100644 --- a/ortools/util/sorted_interval_list.cc +++ b/ortools/util/sorted_interval_list.cc @@ -139,9 +139,12 @@ bool Domain::IsEmpty() const { return intervals_.empty(); } int64 Domain::Size() const { int64 size = 0; for (const ClosedInterval interval : intervals_) { - size += operations_research::CapAdd( - 1, operations_research::CapSub(interval.end, interval.start)); + size = operations_research::CapAdd( + size, operations_research::CapSub(interval.end, interval.start)); } + // Because the intervals are closed on both side above, with miss 1 per + // interval. + size = operations_research::CapAdd(size, intervals_.size()); return size; }