Google OR-Tools v9.15
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
lb_tree_search.cc
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
15
16#include <algorithm>
17#include <cstdint>
18#include <functional>
19#include <limits>
20#include <memory>
21#include <string>
22#include <utility>
23#include <vector>
24
25#include "absl/cleanup/cleanup.h"
26#include "absl/log/check.h"
27#include "absl/log/log.h"
28#include "absl/log/vlog_is_on.h"
29#include "absl/strings/str_cat.h"
30#include "absl/types/span.h"
34#include "ortools/sat/integer.h"
39#include "ortools/sat/model.h"
46#include "ortools/sat/util.h"
49
50namespace operations_research {
51namespace sat {
52
54 : name_(model->Name()),
55 time_limit_(model->GetOrCreate<TimeLimit>()),
56 random_(model->GetOrCreate<ModelRandomGenerator>()),
57 sat_solver_(model->GetOrCreate<SatSolver>()),
58 integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
59 trail_(model->GetOrCreate<Trail>()),
60 assignment_(trail_->Assignment()),
61 integer_trail_(model->GetOrCreate<IntegerTrail>()),
62 watcher_(model->GetOrCreate<GenericLiteralWatcher>()),
63 shared_response_(model->GetOrCreate<SharedResponseManager>()),
64 pseudo_costs_(model->GetOrCreate<PseudoCosts>()),
65 sat_decision_(model->GetOrCreate<SatDecisionPolicy>()),
66 search_helper_(model->GetOrCreate<IntegerSearchHelper>()),
67 parameters_(*model->GetOrCreate<SatParameters>()) {
68 // We should create this class only in the presence of an objective.
69 //
70 // TODO(user): Starts with an initial variable score for all variable in
71 // the objective at their minimum value? this should emulate the first step of
72 // the core approach and gives a similar bound.
73 const ObjectiveDefinition* objective = model->Get<ObjectiveDefinition>();
74 DCHECK(objective != nullptr);
75 objective_var_ = objective->objective_var;
76
77 // Identify an LP with the same objective variable.
78 //
79 // TODO(user): if we have many independent LP, this will find nothing.
82 if (lp->ObjectiveVariable() == objective_var_) {
83 lp_constraint_ = lp;
84 }
85 }
86
87 // We use the normal SAT search but we will bump the variable activity
88 // slightly differently. In addition to the conflicts, we also bump it each
89 // time the objective lower bound increase in a sub-node.
90 std::vector<std::function<BooleanOrIntegerLiteral()>> heuristics;
91 if (SaveLpBasisOption()) {
92 heuristics.emplace_back(LpPseudoCostHeuristic(model));
93 }
94 heuristics.emplace_back(SatSolverHeuristic(model));
95 heuristics.emplace_back(MostFractionalHeuristic(model));
96 heuristics.emplace_back(IntegerValueSelectionHeuristic(
97 model->GetOrCreate<SearchHeuristics>()->fixed_search, model));
98 search_heuristic_ = SequentialSearch(std::move(heuristics));
99}
100
101bool LbTreeSearch::NodeHasBasis(const Node& node) const {
102 return !node.basis.IsEmpty();
103}
104
105bool LbTreeSearch::NodeHasUpToDateBasis(const Node& node) const {
106 if (node.basis.IsEmpty()) return false;
107
108 // TODO(user): Do something smarter. We can at least reuse the variable
109 // statuses maybe?
110 if (node.basis_timestamp != lp_constraint_->num_lp_changes()) {
111 return false;
112 }
113 return true;
114}
115
116void LbTreeSearch::EnableLpAndLoadBestBasis() {
117 DCHECK(lp_constraint_ != nullptr);
118 lp_constraint_->EnablePropagation(true);
119
120 const int level = trail_->CurrentDecisionLevel();
121 if (current_branch_.empty()) return;
122
123 NodeIndex n = current_branch_[0]; // Root.
124 int basis_level = -1;
125 NodeIndex last_node_with_basis(-1);
126 for (int i = 0; i < level; ++i) {
127 if (n >= nodes_.size()) break;
128 if (NodeHasBasis(nodes_[n])) {
129 basis_level = i;
130 last_node_with_basis = n;
131 }
132 const Literal decision = trail_->Decisions()[i].literal;
133 if (nodes_[n].literal_index == decision.Index()) {
134 n = nodes_[n].true_child;
135 } else {
136 DCHECK_EQ(nodes_[n].literal_index, decision.NegatedIndex());
137 n = nodes_[n].false_child;
138 }
139 }
140 if (n < nodes_.size()) {
141 if (NodeHasBasis(nodes_[n])) {
142 basis_level = level;
143 last_node_with_basis = n;
144 }
145 }
146
147 if (last_node_with_basis == -1) {
148 VLOG(1) << "no basis?";
149 return;
150 }
151 VLOG(1) << "load " << basis_level << " / " << level;
152
153 if (!NodeHasUpToDateBasis(nodes_[last_node_with_basis])) {
154 // The basis is no longer up to date, for now we do not load it.
155 // TODO(user): try to do something about it.
156 VLOG(1) << "Skipping potentially bad basis.";
157 return;
158 }
159
160 lp_constraint_->LoadBasisState(nodes_[last_node_with_basis].basis);
161}
162
163void LbTreeSearch::SaveLpBasisInto(Node& node) {
164 node.basis_timestamp = lp_constraint_->num_lp_changes();
165 node.basis = lp_constraint_->GetBasisState();
166}
167
168void LbTreeSearch::UpdateParentObjective(int level) {
169 DCHECK_GE(level, 0);
170 DCHECK_LT(level, current_branch_.size());
171 if (level == 0) return;
172 const NodeIndex parent_index = current_branch_[level - 1];
173 Node& parent = nodes_[parent_index];
174 const NodeIndex child_index = current_branch_[level];
175 const Node& child = nodes_[child_index];
176 if (parent.true_child == child_index) {
177 parent.UpdateTrueObjective(child.MinObjective());
178 } else {
179 DCHECK_EQ(parent.false_child, child_index);
180 parent.UpdateFalseObjective(child.MinObjective());
181 }
182}
183
184void LbTreeSearch::UpdateObjectiveFromParent(int level) {
185 DCHECK_GE(level, 0);
186 DCHECK_LT(level, current_branch_.size());
187 if (level == 0) return;
188 const NodeIndex parent_index = current_branch_[level - 1];
189 const Node& parent = nodes_[parent_index];
190 DCHECK_GE(parent.MinObjective(), current_objective_lb_);
191 const NodeIndex child_index = current_branch_[level];
192 Node& child = nodes_[child_index];
193 if (parent.true_child == child_index) {
194 child.UpdateObjective(parent.true_objective);
195 } else {
196 DCHECK_EQ(parent.false_child, child_index);
197 child.UpdateObjective(parent.false_objective);
198 }
199}
200
201std::string LbTreeSearch::NodeDebugString(NodeIndex n) const {
202 const IntegerValue root_lb = current_objective_lb_;
203 const auto shifted_lb = [root_lb](IntegerValue lb) {
204 return std::max<int64_t>(0, (lb - root_lb).value());
205 };
206
207 std::string s;
208 absl::StrAppend(&s, "#", n.value());
209
210 const Node& node = nodes_[n];
211 std::string true_letter = "t";
212 std::string false_letter = "f";
213 if (node.literal_index != kNoLiteralIndex && !node.is_deleted) {
214 const Literal decision = node.Decision();
215 if (assignment_.LiteralIsTrue(decision)) {
216 true_letter = "T";
217 }
218 if (assignment_.LiteralIsFalse(decision)) {
219 false_letter = "F";
220 }
221 }
222
223 if (node.true_child < nodes_.size()) {
224 absl::StrAppend(&s, " [", true_letter, ":#", node.true_child.value(), " ",
225 shifted_lb(node.true_objective), "]");
226 } else {
227 absl::StrAppend(&s, " [", true_letter, ":## ",
228 shifted_lb(node.true_objective), "]");
229 }
230 if (node.false_child < nodes_.size()) {
231 absl::StrAppend(&s, " [", false_letter, ":#", node.false_child.value(), " ",
232 shifted_lb(node.false_objective), "]");
233 } else {
234 absl::StrAppend(&s, " [", false_letter, ":## ",
235 shifted_lb(node.false_objective), "]");
236 }
237
238 if (node.is_deleted) {
239 absl::StrAppend(&s, " <D>");
240 }
241 if (NodeHasBasis(node)) {
242 absl::StrAppend(&s, " <B>");
243 }
244
245 return s;
246}
247
248void LbTreeSearch::DebugDisplayTree(NodeIndex root) const {
249 int num_nodes = 0;
250
251 util_intops::StrongVector<NodeIndex, int> level(nodes_.size(), 0);
252 std::vector<NodeIndex> to_explore = {root};
253 while (!to_explore.empty()) {
254 NodeIndex n = to_explore.back();
255 to_explore.pop_back();
256
257 ++num_nodes;
258 const Node& node = nodes_[n];
259
260 std::string s(level[n], ' ');
261 absl::StrAppend(&s, NodeDebugString(n));
262 LOG(INFO) << s;
263
264 if (node.true_child < nodes_.size()) {
265 to_explore.push_back(node.true_child);
266 level[node.true_child] = level[n] + 1;
267 }
268 if (node.false_child < nodes_.size()) {
269 to_explore.push_back(node.false_child);
270 level[node.false_child] = level[n] + 1;
271 }
272 }
273 LOG(INFO) << "num_nodes: " << num_nodes;
274}
275
276// Here we forgot the whole search tree and restart.
277//
278// The idea is that the heuristic has now more information so it will likely
279// take better decision which will result in a smaller overall tree.
280bool LbTreeSearch::FullRestart() {
281 ++num_full_restarts_;
282 num_decisions_taken_at_last_restart_ = num_decisions_taken_;
283 num_nodes_in_tree_ = 0;
284 nodes_.clear();
285 current_branch_.clear();
286 return sat_solver_->ResetToLevelZero();
287}
288
289void LbTreeSearch::MarkAsDeletedNodeAndUnreachableSubtree(Node& node) {
290 --num_nodes_in_tree_;
291 DCHECK(!node.is_deleted);
292 node.is_deleted = true;
293 DCHECK_NE(node.literal_index, kNoLiteralIndex);
294 if (assignment_.LiteralIsTrue(Literal(node.literal_index))) {
295 MarkSubtreeAsDeleted(node.false_child);
296 } else {
297 MarkSubtreeAsDeleted(node.true_child);
298 }
299}
300
301void LbTreeSearch::MarkBranchAsInfeasible(Node& node, bool true_branch) {
302 if (true_branch) {
303 node.UpdateTrueObjective(kMaxIntegerValue);
304 MarkSubtreeAsDeleted(node.true_child);
305 node.true_child = NodeIndex(std::numeric_limits<int32_t>::max());
306 } else {
307 node.UpdateFalseObjective(kMaxIntegerValue);
308 MarkSubtreeAsDeleted(node.false_child);
309 node.false_child = NodeIndex(std::numeric_limits<int32_t>::max());
310 }
311}
312
313void LbTreeSearch::MarkSubtreeAsDeleted(NodeIndex root) {
314 std::vector<NodeIndex> to_delete{root};
315 for (int i = 0; i < to_delete.size(); ++i) {
316 const NodeIndex n = to_delete[i];
317 if (n >= nodes_.size()) continue;
318 if (nodes_[n].is_deleted) continue;
319
320 --num_nodes_in_tree_;
321 nodes_[n].is_deleted = true;
322
323 to_delete.push_back(nodes_[n].true_child);
324 to_delete.push_back(nodes_[n].false_child);
325 }
326}
327
328std::string LbTreeSearch::SmallProgressString() const {
329 return absl::StrCat(
330 "nodes=", num_nodes_in_tree_, "/", nodes_.size(),
331 " rc=", num_rc_detected_, " decisions=", num_decisions_taken_,
332 " @root=", num_back_to_root_node_, " restarts=", num_full_restarts_,
333 " lp_iters=[", FormatCounter(num_lp_iters_at_level_zero_), ", ",
334 FormatCounter(num_lp_iters_save_basis_), ", ",
335 FormatCounter(num_lp_iters_first_branch_), ", ",
336 FormatCounter(num_lp_iters_dive_), "]");
337}
338
339std::function<void()> LbTreeSearch::UpdateLpIters(int64_t* counter) {
340 if (lp_constraint_ == nullptr) return []() {};
341 const int64_t old_num = lp_constraint_->total_num_simplex_iterations();
342 return [old_num, counter, this]() {
343 const int64_t new_num = lp_constraint_->total_num_simplex_iterations();
344 *counter += new_num - old_num;
345 };
346}
347
348bool LbTreeSearch::LevelZeroLogic() {
349 ++num_back_to_root_node_;
350 num_decisions_taken_at_last_level_zero_ = num_decisions_taken_;
351
352 // Always run the LP when we are back at level zero.
353 if (SaveLpBasisOption() && !current_branch_.empty()) {
354 const auto cleanup =
355 absl::MakeCleanup(UpdateLpIters(&num_lp_iters_at_level_zero_));
356 EnableLpAndLoadBestBasis();
357 if (!sat_solver_->FinishPropagation()) {
358 return false;
359 }
360 SaveLpBasisInto(nodes_[current_branch_[0]]);
361 lp_constraint_->EnablePropagation(false);
362 }
363
364 // Import the objective upper-bound.
365 // We do that manually because we disabled objective import to not "pollute"
366 // the objective lower_bound and still have local reason for objective
367 // improvement.
368 {
369 const IntegerValue ub = shared_response_->GetInnerObjectiveUpperBound();
370 if (integer_trail_->UpperBound(objective_var_) > ub) {
371 if (!integer_trail_->Enqueue(
372 IntegerLiteral::LowerOrEqual(objective_var_, ub), {}, {})) {
373 sat_solver_->NotifyThatModelIsUnsat();
374 return false;
375 }
376 if (!sat_solver_->FinishPropagation()) {
377 return false;
378 }
379 }
380 }
381
382 // If the search has not just been restarted (in which case nodes_ would be
383 // empty), and if we are at level zero (either naturally, or if the
384 // backtrack level was set to zero in the above code), let's run a different
385 // heuristic to decide whether to restart the search from scratch or not.
386 //
387 // We ignore small search trees.
388 if (num_nodes_in_tree_ > 50) {
389 // Let's count how many nodes have worse objective bounds than the best
390 // known external objective lower bound.
391 const IntegerValue latest_lb =
392 shared_response_->GetInnerObjectiveLowerBound();
393 int num_nodes = 0;
394 int num_nodes_with_lower_objective = 0;
395 for (const Node& node : nodes_) {
396 if (node.is_deleted) continue;
397 ++num_nodes;
398 if (node.MinObjective() < latest_lb) num_nodes_with_lower_objective++;
399 }
400 DCHECK_EQ(num_nodes_in_tree_, num_nodes);
401 if (num_nodes_with_lower_objective * 2 > num_nodes) {
402 VLOG(2) << "lb_tree_search restart nodes: "
403 << num_nodes_with_lower_objective << "/" << num_nodes << " : "
404 << 100.0 * num_nodes_with_lower_objective / num_nodes << "%"
405 << ", decisions:" << num_decisions_taken_;
406 if (!FullRestart()) return false;
407 }
408 }
409
410 return true;
411}
412
414 const std::function<void()>& feasible_solution_observer) {
415 if (!sat_solver_->ResetToLevelZero()) {
416 return sat_solver_->UnsatStatus();
417 }
418
419 // We currently restart the search tree from scratch from time to times:
420 // - Initially, every kNumDecisionsBeforeInitialRestarts, for at most
421 // kMaxNumInitialRestarts times.
422 // - Every time we backtrack to level zero, we count how many nodes are worse
423 // than the best known objective lower bound. If this is true for more than
424 // half of the existing nodes, we restart and clear all nodes. If if this
425 // happens during the initial restarts phase, it reset the above counter and
426 // uses 1 of the available initial restarts.
427 //
428 // This has 2 advantages:
429 // - It allows our "pseudo-cost" to kick in and experimentally result in
430 // smaller trees down the road.
431 // - It removes large inefficient search trees.
432 //
433 // TODO(user): a strong branching initial start, or allowing a few decision
434 // per nodes might be a better approach.
435 //
436 // TODO(user): It would also be cool to exploit the reason for the LB increase
437 // even more.
438 const int kMaxNumInitialRestarts = 10;
439 const int64_t kNumDecisionsBeforeInitialRestarts = 1000;
440
441 // If some branches already have a good lower bound, no need to call the LP
442 // on those.
443 watcher_->SetStopPropagationCallback([this] {
444 return integer_trail_->LowerBound(objective_var_) > current_objective_lb_;
445 });
446
447 while (!time_limit_->LimitReached() && !shared_response_->ProblemIsSolved()) {
448 VLOG(2) << "LOOP " << sat_solver_->CurrentDecisionLevel();
449
450 // Each time we are back here, we bump the activities of the variable that
451 // are part of the objective lower bound reason.
452 //
453 // Note that this is why we prefer not to increase the lower zero lower
454 // bound of objective_var_ with the tree root lower bound, so we can exploit
455 // more reasons.
456 //
457 // TODO(user): This is slightly different than bumping each time we
458 // push a decision that result in an LB increase. This is also called on
459 // backjump for instance.
460 const IntegerValue obj_diff =
461 integer_trail_->LowerBound(objective_var_) -
462 integer_trail_->LevelZeroLowerBound(objective_var_);
463 if (obj_diff > 0) {
464 std::vector<Literal> reason =
465 integer_trail_->ReasonFor(IntegerLiteral::GreaterOrEqual(
466 objective_var_, integer_trail_->LowerBound(objective_var_)));
467
468 // TODO(user): We also need to update pseudo cost on conflict.
469 pseudo_costs_->UpdateBoolPseudoCosts(reason, obj_diff);
470 sat_decision_->BumpVariableActivities(reason);
471 sat_decision_->UpdateVariableActivityIncrement();
472 }
473
474 // Propagate upward in the tree the new objective lb.
475 if (!current_branch_.empty()) {
476 // Our branch is always greater or equal to the level.
477 // We increase the objective_lb of the current node if needed.
478 {
479 const int current_level = sat_solver_->CurrentDecisionLevel();
480 const IntegerValue current_objective_lb =
481 integer_trail_->LowerBound(objective_var_);
482 if (DEBUG_MODE) {
483 CHECK_LE(current_level, current_branch_.size());
484 for (int i = 0; i < current_level; ++i) {
485 CHECK(!nodes_[current_branch_[i]].is_deleted);
486 CHECK(assignment_.LiteralIsAssigned(
487 nodes_[current_branch_[i]].Decision()));
488 }
489 }
490 if (current_level < current_branch_.size()) {
491 nodes_[current_branch_[current_level]].UpdateObjective(
492 current_objective_lb);
493 }
494
495 // Minor optim: sometimes, because of the LP and cuts, the reason for
496 // objective_var_ only contains lower level literals, so we can exploit
497 // that.
498 //
499 // TODO(user): No point checking that if the objective lb wasn't
500 // assigned at this level.
501 //
502 // TODO(user): Exploit the reasons further.
503 if (current_objective_lb >
504 integer_trail_->LevelZeroLowerBound(objective_var_)) {
505 const std::vector<Literal> reason =
506 integer_trail_->ReasonFor(IntegerLiteral::GreaterOrEqual(
507 objective_var_, current_objective_lb));
508 int max_level = 0;
509 for (const Literal l : reason) {
510 max_level = std::max<int>(
511 max_level,
512 sat_solver_->LiteralTrail().Info(l.Variable()).level);
513 }
514 if (max_level < current_level) {
515 nodes_[current_branch_[max_level]].UpdateObjective(
516 current_objective_lb);
517 }
518 }
519 }
520
521 // Propagate upward any new bounds.
522 for (int level = current_branch_.size(); --level > 0;) {
523 UpdateParentObjective(level);
524 }
525 }
526
527 if (SaveLpBasisOption()) {
528 // We disable LP automatic propagation and only enable it:
529 // - at root node
530 // - when we go to a new branch.
531 lp_constraint_->EnablePropagation(false);
532 }
533 if (sat_solver_->ModelIsUnsat()) return sat_solver_->UnsatStatus();
534
535 // This will import other workers bound if we are back to level zero.
536 // It might also decide to restart.
537 if (!search_helper_->BeforeTakingDecision()) {
538 return sat_solver_->UnsatStatus();
539 }
540
541 // This is the current bound we try to improve. We cache it here to avoid
542 // getting the lock many times and it is also easier to follow the code if
543 // this is assumed constant for one iteration.
544 current_objective_lb_ = shared_response_->GetInnerObjectiveLowerBound();
545 if (!current_branch_.empty()) {
546 nodes_[current_branch_[0]].UpdateObjective(current_objective_lb_);
547 for (int i = 1; i < current_branch_.size(); ++i) {
548 UpdateObjectiveFromParent(i);
549 }
550
551 // If the root lb increased, update global shared objective lb.
552 const IntegerValue bound = nodes_[current_branch_[0]].MinObjective();
553 if (bound > current_objective_lb_) {
554 shared_response_->UpdateInnerObjectiveBounds(
555 absl::StrCat(name_, " (", SmallProgressString(), ") "), bound,
556 integer_trail_->LevelZeroUpperBound(objective_var_));
557 current_objective_lb_ = bound;
558 if (VLOG_IS_ON(3)) DebugDisplayTree(current_branch_[0]);
559 }
560 }
561
562 // Forget the whole tree and restart.
563 // We will do it periodically at the beginning of the search each time we
564 // cross the kNumDecisionsBeforeInitialRestarts decision since the last
565 // restart. This will happen at most kMaxNumInitialRestarts times.
566 if (num_decisions_taken_ >= num_decisions_taken_at_last_restart_ +
567 kNumDecisionsBeforeInitialRestarts &&
568 num_full_restarts_ < kMaxNumInitialRestarts) {
569 VLOG(2) << "lb_tree_search (initial_restart " << SmallProgressString()
570 << ")";
571 if (!FullRestart()) return sat_solver_->UnsatStatus();
572 continue;
573 }
574
575 // Periodic backtrack to level zero so we can import bounds.
576 if (num_decisions_taken_ >=
577 num_decisions_taken_at_last_level_zero_ + 10000) {
578 if (!sat_solver_->ResetToLevelZero()) {
579 return sat_solver_->UnsatStatus();
580 }
581 }
582
583 // Backtrack if needed.
584 //
585 // Our algorithm stop exploring a branch as soon as its objective lower
586 // bound is greater than the root lower bound. We then backtrack to the
587 // first node in the branch that is not yet closed under this bound.
588 //
589 // TODO(user): If we remember how far we can backjump for both true/false
590 // branch, we could be more efficient.
591 while (current_branch_.size() > sat_solver_->CurrentDecisionLevel() + 1 ||
592 (current_branch_.size() > 1 &&
593 nodes_[current_branch_.back()].MinObjective() >
594 current_objective_lb_)) {
595 current_branch_.pop_back();
596 }
597
598 // Backtrack the solver to be in sync with current_branch_.
599 {
600 const int backtrack_level =
601 std::max(0, static_cast<int>(current_branch_.size()) - 1);
602 sat_solver_->Backtrack(backtrack_level);
603 if (!sat_solver_->FinishPropagation()) {
604 return sat_solver_->UnsatStatus();
605 }
606 if (sat_solver_->CurrentDecisionLevel() < backtrack_level) {
607 continue;
608 }
609 }
610
611 if (sat_solver_->CurrentDecisionLevel() == 0) {
612 if (!LevelZeroLogic()) {
613 return sat_solver_->UnsatStatus();
614 }
615 }
616
617 // Dive: Follow the branch with lowest objective.
618 // Note that we do not creates new nodes here.
619 //
620 // TODO(user): If we have new information and our current objective bound
621 // is higher than any bound in a whole subtree, we might want to just
622 // restart this subtree exploration?
623 while (true) {
624 const int size = current_branch_.size();
625 const int level = sat_solver_->CurrentDecisionLevel();
626
627 // Invariant are tricky:
628 // current_branch_ contains one entry per decision taken + the last one
629 // which we are about to take. If we don't have the last entry, it means
630 // we are about to take a new decision.
631 DCHECK(size == level || size == level + 1);
632 if (size == level) break; // Take new decision.
633
634 const NodeIndex node_index = current_branch_[level];
635 Node& node = nodes_[node_index];
636 DCHECK_GT(node.true_child, node_index);
637 DCHECK_GT(node.false_child, node_index);
638
639 // If the bound of this node is high, restart the main loop..
640 node.UpdateObjective(std::max(
641 current_objective_lb_, integer_trail_->LowerBound(objective_var_)));
642 if (node.MinObjective() > current_objective_lb_) break;
643 DCHECK_EQ(node.MinObjective(), current_objective_lb_) << level;
644
645 // This will be set to the next node index.
646 NodeIndex n;
647 DCHECK(!node.is_deleted);
648 const Literal node_literal = node.Decision();
649
650 // If the variable is already fixed, we bypass the node and connect
651 // its parent directly to the relevant child.
652 if (assignment_.LiteralIsAssigned(node_literal)) {
653 IntegerValue new_lb;
654 if (assignment_.LiteralIsTrue(node_literal)) {
655 n = node.true_child;
656 new_lb = node.true_objective;
657 } else {
658 n = node.false_child;
659 new_lb = node.false_objective;
660 }
661 MarkAsDeletedNodeAndUnreachableSubtree(node);
662
663 // We jump directly to the subnode.
664 // Else we will change the root.
665 current_branch_.pop_back();
666 if (!current_branch_.empty()) {
667 const NodeIndex parent = current_branch_.back();
668 DCHECK(!nodes_[parent].is_deleted);
669 const Literal parent_literal = nodes_[parent].Decision();
670 if (assignment_.LiteralIsTrue(parent_literal)) {
671 nodes_[parent].true_child = n;
672 nodes_[parent].UpdateTrueObjective(new_lb);
673 } else {
674 DCHECK(assignment_.LiteralIsFalse(parent_literal));
675 nodes_[parent].false_child = n;
676 nodes_[parent].UpdateFalseObjective(new_lb);
677 }
678 if (new_lb > current_objective_lb_) {
679 // This is probably not needed.
680 if (n < nodes_.size() && !nodes_[n].IsLeaf()) {
681 current_branch_.push_back(n);
682 nodes_[n].UpdateObjective(new_lb);
683 }
684 break;
685 }
686 } else {
687 if (n >= nodes_.size()) {
688 // We never explored the other branch, so we can just clear all
689 // nodes.
690 num_nodes_in_tree_ = 0;
691 nodes_.clear();
692 } else if (nodes_[n].IsLeaf()) {
693 // Keep the saved basis.
694 num_nodes_in_tree_ = 1;
695 Node root = std::move(nodes_[n]);
696 nodes_.clear();
697 nodes_.push_back(std::move(root));
698 n = NodeIndex(0);
699 } else {
700 // We always make sure the root is at zero.
701 // The root is no longer at zero, that might cause issue.
702 // Cleanup.
703 }
704 }
705 } else {
706 // See if we have better bounds using the current LP state.
707 ExploitReducedCosts(current_branch_[level]);
708 if (node.MinObjective() > current_objective_lb_) break;
709
710 // If both lower bound are the same, we pick the literal branch. We do
711 // that because this is the polarity that was chosen by the SAT
712 // heuristic in the first place. We tried random, it doesn't seems to
713 // work as well.
714 num_decisions_taken_++;
715 const bool choose_true = node.true_objective <= node.false_objective;
716 Literal next_decision;
717 if (choose_true) {
718 n = node.true_child;
719 next_decision = node_literal;
720 } else {
721 n = node.false_child;
722 next_decision = node_literal.Negated();
723 }
724
725 // If we are taking this branch for the first time, we enable the LP and
726 // make sure we solve it before taking the decision. This allow to have
727 // proper pseudo-costs, and also be incremental for the decision we are
728 // about to take.
729 //
730 // We also enable the LP if we have no basis info for this node.
731 if (SaveLpBasisOption() &&
732 (n >= nodes_.size() || !NodeHasBasis(node))) {
733 const auto cleanup =
734 absl::MakeCleanup(UpdateLpIters(&num_lp_iters_save_basis_));
735
736 VLOG(1) << "~~~~";
737 EnableLpAndLoadBestBasis();
738 const int level = sat_solver_->CurrentDecisionLevel();
739 if (!sat_solver_->FinishPropagation()) {
740 return sat_solver_->UnsatStatus();
741 }
742 if (sat_solver_->CurrentDecisionLevel() < level) {
743 node.UpdateObjective(kMaxIntegerValue);
744 break;
745 }
746
747 // The decision might have become assigned, in which case we loop.
748 if (assignment_.LiteralIsAssigned(next_decision)) {
749 continue;
750 }
751
752 SaveLpBasisInto(node);
753
754 // If we are not at the end, disable the LP propagation.
755 if (n < nodes_.size()) {
756 lp_constraint_->EnablePropagation(false);
757 }
758 }
759
760 // Take the decision.
761 const auto cleanup =
762 absl::MakeCleanup(UpdateLpIters(&num_lp_iters_first_branch_));
763 DCHECK(!assignment_.LiteralIsAssigned(next_decision));
764 search_helper_->TakeDecision(next_decision);
765
766 // Conflict?
767 if (current_branch_.size() != sat_solver_->CurrentDecisionLevel()) {
768 MarkBranchAsInfeasible(node, choose_true);
769 break;
770 }
771
772 // Update the proper field and abort the dive if we crossed the
773 // threshold.
774 const IntegerValue lb = integer_trail_->LowerBound(objective_var_);
775 if (choose_true) {
776 node.UpdateTrueObjective(lb);
777 } else {
778 node.UpdateFalseObjective(lb);
779 }
780
781 if (n < nodes_.size()) {
782 nodes_[n].UpdateObjective(lb);
783 } else if (SaveLpBasisOption()) {
784 SaveLpBasisInto(nodes_[CreateNewEmptyNodeIfNeeded()]);
785 }
786
787 if (lb > current_objective_lb_) break;
788 }
789
790 if (VLOG_IS_ON(1)) {
791 shared_response_->LogMessageWithThrottling(
792 "TreeS", absl::StrCat(" (", SmallProgressString(), ")"));
793 }
794
795 if (n < nodes_.size() && !nodes_[n].IsLeaf()) {
796 current_branch_.push_back(n);
797 } else {
798 break;
799 }
800 }
801
802 // If a conflict occurred, we will backtrack.
803 if (current_branch_.size() != sat_solver_->CurrentDecisionLevel()) {
804 continue;
805 }
806
807 // TODO(user): The code is hard to follow. Fix and merge that with test
808 // below.
809 if (!current_branch_.empty()) {
810 const Node& final_node = nodes_[current_branch_.back()];
811 if (assignment_.LiteralIsTrue(final_node.Decision())) {
812 if (final_node.true_objective > current_objective_lb_) {
813 continue;
814 }
815 } else {
816 DCHECK(assignment_.LiteralIsFalse(final_node.Decision()));
817 if (final_node.false_objective > current_objective_lb_) {
818 continue;
819 }
820 }
821 }
822
823 // This test allow to not take a decision when the branch is already closed
824 // (i.e. the true branch or false branch lb is high enough). Adding it
825 // basically changes if we take the decision later when we explore the
826 // branch or right now.
827 //
828 // I feel taking it later is better. It also avoid creating unneeded nodes.
829 // It does change the behavior on a few problem though. For instance on
830 // irp.mps.gz, the search works better without this, whatever the random
831 // seed. Not sure why, maybe it creates more diversity?
832 //
833 // Another difference is that if the search is done and we have a feasible
834 // solution, we will not report it because of this test (except if we are
835 // at the optimal).
836 if (integer_trail_->LowerBound(objective_var_) > current_objective_lb_) {
837 continue;
838 }
839
840 const auto cleanup = absl::MakeCleanup(UpdateLpIters(&num_lp_iters_dive_));
841
842 if (current_branch_.empty()) {
843 VLOG(2) << "DIVE from empty tree";
844 } else {
845 VLOG(2) << "DIVE from " << NodeDebugString(current_branch_.back());
846 }
847
848 if (SaveLpBasisOption() && !lp_constraint_->PropagationIsEnabled()) {
849 // This reuse or create a node to store the basis.
850 const NodeIndex index = CreateNewEmptyNodeIfNeeded();
851
852 EnableLpAndLoadBestBasis();
853 const int level = sat_solver_->CurrentDecisionLevel();
854 if (!sat_solver_->FinishPropagation()) {
855 return sat_solver_->UnsatStatus();
856 }
857
858 // Loop on backtrack or bound improvement.
859 if (sat_solver_->CurrentDecisionLevel() < level) {
860 Node& node = nodes_[index];
861 node.UpdateObjective(kMaxIntegerValue);
862 continue;
863 }
864
865 SaveLpBasisInto(nodes_[index]);
866
867 const IntegerValue obj_lb = integer_trail_->LowerBound(objective_var_);
868 if (obj_lb > current_objective_lb_) {
869 nodes_[index].UpdateObjective(obj_lb);
870 if (!current_branch_.empty()) {
871 Node& parent_node = nodes_[current_branch_.back()];
872 const Literal node_literal = parent_node.Decision();
873 DCHECK(assignment_.LiteralIsAssigned(node_literal));
874 if (assignment_.LiteralIsTrue(node_literal)) {
875 parent_node.UpdateTrueObjective(obj_lb);
876 } else {
877 parent_node.UpdateFalseObjective(obj_lb);
878 }
879 }
880 continue;
881 }
882 }
883
884 // Invariant: The current branch is fully assigned, and the solver is in
885 // sync. And we are not on a "bad" path.
886 const int base_level = sat_solver_->CurrentDecisionLevel();
887 if (DEBUG_MODE) {
888 CHECK_EQ(base_level, current_branch_.size());
889 for (const NodeIndex index : current_branch_) {
890 CHECK(!nodes_[index].is_deleted);
891 const Literal decision = nodes_[index].Decision();
892 if (assignment_.LiteralIsTrue(decision)) {
893 CHECK_EQ(nodes_[index].true_objective, current_objective_lb_);
894 } else {
895 CHECK(assignment_.LiteralIsFalse(decision));
896 CHECK_EQ(nodes_[index].false_objective, current_objective_lb_);
897 }
898 }
899 }
900
901 // We are about to take a new decision, what we will do is dive until
902 // the objective lower bound increase. we will then create a bunch of new
903 // nodes in the tree.
904 //
905 // By analyzing the reason for the increase, we can create less nodes than
906 // if we just followed the initial heuristic.
907 //
908 // TODO(user): In multithread, this change the behavior a lot since we
909 // dive until we beat the best shared bound. Maybe we shouldn't do that.
910 while (true) {
911 // TODO(user): We sometimes branch on the objective variable, this should
912 // probably be avoided.
913 if (sat_solver_->ModelIsUnsat()) return sat_solver_->UnsatStatus();
914 LiteralIndex decision;
915 if (!search_helper_->GetDecision(search_heuristic_, &decision)) {
916 continue;
917 }
918
919 // No new decision: search done.
920 if (time_limit_->LimitReached()) return SatSolver::LIMIT_REACHED;
921 if (decision == kNoLiteralIndex) {
922 feasible_solution_observer();
923 break;
924 }
925
926 num_decisions_taken_++;
927 if (!search_helper_->TakeDecision(Literal(decision))) {
928 return sat_solver_->UnsatStatus();
929 }
930 if (trail_->CurrentDecisionLevel() < base_level) {
931 // TODO(user): it would be nice to mark some node as infeasible if
932 // this is the case. However this could happen after many decision and
933 // we realize with the lp that one of them should have been fixed
934 // earlier, without any infeasibility in the current branch.
935 break;
936 }
937 if (integer_trail_->LowerBound(objective_var_) > current_objective_lb_) {
938 break;
939 }
940 }
941
942 if (trail_->CurrentDecisionLevel() <= base_level) {
943 continue;
944 }
945
946 // Analyse the reason for objective increase. Deduce a set of new nodes to
947 // append to the tree.
948 //
949 // TODO(user): Try to minimize the number of decisions?
950 const std::vector<Literal> reason =
951 integer_trail_->ReasonFor(IntegerLiteral::GreaterOrEqual(
952 objective_var_, integer_trail_->LowerBound(objective_var_)));
953 std::vector<Literal> decisions = ExtractDecisions(base_level, reason);
954
955 // Bump activities.
956 sat_decision_->BumpVariableActivities(reason);
957 sat_decision_->BumpVariableActivities(decisions);
958 sat_decision_->UpdateVariableActivityIncrement();
959
960 // Create one node per new decisions.
961 DCHECK_EQ(current_branch_.size(), base_level);
962 for (const Literal d : decisions) {
963 AppendNewNodeToCurrentBranch(d);
964 }
965
966 // TODO(user): We should probably save the basis in more cases.
967 if (SaveLpBasisOption() && decisions.size() == 1) {
968 SaveLpBasisInto(nodes_[CreateNewEmptyNodeIfNeeded()]);
969 }
970
971 // Update the objective of the last node in the branch since we just
972 // improved that.
973 if (!current_branch_.empty()) {
974 Node& n = nodes_[current_branch_.back()];
975 if (assignment_.LiteralIsTrue(n.Decision())) {
976 n.UpdateTrueObjective(integer_trail_->LowerBound(objective_var_));
977 } else {
978 n.UpdateFalseObjective(integer_trail_->LowerBound(objective_var_));
979 }
980 }
981
982 // Reset the solver to a correct state since we have a subset of the
983 // current propagation. We backtrack as little as possible.
984 //
985 // The decision level is the number of decision taken.
986 // Decision()[level] is the decision at that level.
987 int backtrack_level = base_level;
988 DCHECK_LE(current_branch_.size(), trail_->CurrentDecisionLevel());
989 while (backtrack_level < current_branch_.size() &&
990 trail_->Decisions()[backtrack_level].literal.Index() ==
991 nodes_[current_branch_[backtrack_level]].literal_index) {
992 ++backtrack_level;
993 }
994 sat_solver_->BacktrackAndPropagateReimplications(backtrack_level);
995
996 // Update bounds with reduced costs info.
997 //
998 // TODO(user): Uses old optimal constraint that we just potentially
999 // backtracked over?
1000 //
1001 // TODO(user): We could do all at once rather than in O(#decision * #size).
1002 for (int i = backtrack_level; i < current_branch_.size(); ++i) {
1003 ExploitReducedCosts(current_branch_[i]);
1004 }
1005 }
1006
1008}
1009
1010std::vector<Literal> LbTreeSearch::ExtractDecisions(
1011 int base_level, absl::Span<const Literal> conflict) {
1012 std::vector<int> num_per_level(sat_solver_->CurrentDecisionLevel() + 1, 0);
1013 std::vector<bool> is_marked;
1014 for (const Literal l : conflict) {
1015 const AssignmentInfo& info = trail_->Info(l.Variable());
1016 if (info.level <= base_level) continue;
1017 num_per_level[info.level]++;
1018 if (info.trail_index >= is_marked.size()) {
1019 is_marked.resize(info.trail_index + 1);
1020 }
1021 is_marked[info.trail_index] = true;
1022 }
1023
1024 std::vector<Literal> result;
1025 if (is_marked.empty()) return result;
1026 for (int i = is_marked.size() - 1; i >= 0; --i) {
1027 if (!is_marked[i]) continue;
1028
1029 const Literal l = (*trail_)[i];
1030 const AssignmentInfo& info = trail_->Info(l.Variable());
1031 if (info.level <= base_level) break;
1032 if (num_per_level[info.level] == 1) {
1033 result.push_back(l);
1034 continue;
1035 }
1036
1037 // Expand.
1038 num_per_level[info.level]--;
1039 for (const Literal new_l : trail_->Reason(l.Variable())) {
1040 const AssignmentInfo& new_info = trail_->Info(new_l.Variable());
1041 if (new_info.level <= base_level) continue;
1042 if (is_marked[new_info.trail_index]) continue;
1043 is_marked[new_info.trail_index] = true;
1044 num_per_level[new_info.level]++;
1045 }
1046 }
1047
1048 // We prefer to keep the same order.
1049 std::reverse(result.begin(), result.end());
1050 return result;
1051}
1052
1053LbTreeSearch::NodeIndex LbTreeSearch::CreateNewEmptyNodeIfNeeded() {
1054 NodeIndex n(0);
1055 if (current_branch_.empty()) {
1056 if (nodes_.empty()) {
1057 ++num_nodes_in_tree_;
1058 nodes_.emplace_back(current_objective_lb_);
1059 } else {
1060 DCHECK_EQ(nodes_.size(), 1);
1061 }
1062 } else {
1063 const NodeIndex parent = current_branch_.back();
1064 DCHECK(!nodes_[parent].is_deleted);
1065 const Literal parent_literal = nodes_[parent].Decision();
1066 if (assignment_.LiteralIsTrue(parent_literal)) {
1067 if (nodes_[parent].true_child >= nodes_.size()) {
1068 n = NodeIndex(nodes_.size());
1069 ++num_nodes_in_tree_;
1070 nodes_[parent].true_child = NodeIndex(nodes_.size());
1071 nodes_.emplace_back(current_objective_lb_);
1072 } else {
1073 n = nodes_[parent].true_child;
1074 }
1075 nodes_[parent].UpdateTrueObjective(current_objective_lb_);
1076 } else {
1077 DCHECK(assignment_.LiteralIsFalse(parent_literal));
1078 if (nodes_[parent].false_child >= nodes_.size()) {
1079 n = NodeIndex(nodes_.size());
1080 ++num_nodes_in_tree_;
1081 nodes_[parent].false_child = NodeIndex(nodes_.size());
1082 nodes_.emplace_back(current_objective_lb_);
1083 } else {
1084 n = nodes_[parent].false_child;
1085 }
1086 nodes_[parent].UpdateFalseObjective(current_objective_lb_);
1087 }
1088 }
1089 DCHECK(!nodes_[n].is_deleted);
1090 DCHECK_EQ(nodes_[n].literal_index, kNoLiteralIndex);
1091 return n;
1092}
1093
1094void LbTreeSearch::AppendNewNodeToCurrentBranch(Literal decision) {
1095 NodeIndex n(nodes_.size());
1096 if (current_branch_.empty()) {
1097 if (nodes_.empty()) {
1098 ++num_nodes_in_tree_;
1099 nodes_.emplace_back(current_objective_lb_);
1100 } else {
1101 DCHECK_EQ(nodes_.size(), 1);
1102 n = 0;
1103 }
1104 } else {
1105 const NodeIndex parent = current_branch_.back();
1106 DCHECK(!nodes_[parent].is_deleted);
1107 const Literal parent_literal = nodes_[parent].Decision();
1108 if (assignment_.LiteralIsTrue(parent_literal)) {
1109 if (nodes_[parent].true_child < nodes_.size()) {
1110 n = nodes_[parent].true_child;
1111 } else {
1112 ++num_nodes_in_tree_;
1113 nodes_.emplace_back(current_objective_lb_);
1114 nodes_[parent].true_child = n;
1115 }
1116 nodes_[parent].UpdateTrueObjective(current_objective_lb_);
1117 } else {
1118 DCHECK(assignment_.LiteralIsFalse(parent_literal));
1119 if (nodes_[parent].false_child < nodes_.size()) {
1120 n = nodes_[parent].false_child;
1121 } else {
1122 ++num_nodes_in_tree_;
1123 nodes_.emplace_back(current_objective_lb_);
1124 nodes_[parent].false_child = n;
1125 }
1126 nodes_[parent].UpdateFalseObjective(current_objective_lb_);
1127 }
1128 }
1129 DCHECK_LT(n, nodes_.size());
1130 DCHECK_EQ(nodes_[n].literal_index, kNoLiteralIndex) << " issue " << n;
1131 nodes_[n].SetDecision(decision);
1132 nodes_[n].UpdateObjective(current_objective_lb_);
1133 current_branch_.push_back(n);
1134}
1135
1136// Looking at the reduced costs, we can already have a bound for one of the
1137// branch. Increasing the corresponding objective can save some branches,
1138// and also allow for a more incremental LP solving since we do less back
1139// and forth.
1140//
1141// TODO(user): The code to recover that is a bit convoluted. Alternatively
1142// Maybe we should do a "fast" propagation without the LP in each branch.
1143// That will work as long as we keep these optimal LP constraints around
1144// and propagate them.
1145//
1146// TODO(user): Incorporate this in the heuristic so we choose more Booleans
1147// inside these LP explanations?
1148void LbTreeSearch::ExploitReducedCosts(NodeIndex n) {
1149 if (lp_constraint_ == nullptr) return;
1150
1151 // TODO(user): we could consider earlier constraints instead of just
1152 // looking at the last one, but experiments didn't really show a big
1153 // gain.
1154 const auto& cts = lp_constraint_->OptimalConstraints();
1155 if (cts.empty()) return;
1156 const std::unique_ptr<IntegerSumLE128>& rc = cts.back();
1157
1158 // Note that this return literal EQUIVALENT to the node.literal, not just
1159 // implied by it. We need that for correctness.
1160 int num_tests = 0;
1161 Node& node = nodes_[n];
1162 DCHECK(!node.is_deleted);
1163 const Literal node_literal = node.Decision();
1164
1165 // This can happen if we have re-implication and propagation...
1166 if (assignment_.LiteralIsAssigned(node_literal)) return;
1167
1168 for (const IntegerLiteral integer_literal :
1169 integer_encoder_->GetIntegerLiterals(node_literal)) {
1170 // To avoid bad corner case. Not sure it ever triggers.
1171 if (++num_tests > 10) break;
1172
1173 const std::pair<IntegerValue, IntegerValue> bounds =
1174 rc->ConditionalLb(integer_literal, objective_var_);
1175 if (bounds.first > node.false_objective) {
1176 ++num_rc_detected_;
1177 node.UpdateFalseObjective(bounds.first);
1178 }
1179 if (bounds.second > node.true_objective) {
1180 ++num_rc_detected_;
1181 node.UpdateTrueObjective(bounds.second);
1182 }
1183 }
1184}
1185
1186} // namespace sat
1187} // namespace operations_research
SatSolver::Status Search(const std::function< void()> &feasible_solution_observer)
T Get(std::function< T(const Model &)> f) const
Similar to Add() but this is const.
Definition model.h:110
const AssignmentInfo & Info(BooleanVariable var) const
Definition sat_base.h:655
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
const LiteralIndex kNoLiteralIndex(-1)
std::function< BooleanOrIntegerLiteral()> IntegerValueSelectionHeuristic(std::function< BooleanOrIntegerLiteral()> var_selection_heuristic, Model *model)
std::function< BooleanOrIntegerLiteral()> SatSolverHeuristic(Model *model)
std::function< BooleanOrIntegerLiteral()> LpPseudoCostHeuristic(Model *model)
std::function< BooleanOrIntegerLiteral()> MostFractionalHeuristic(Model *model)
std::function< BooleanOrIntegerLiteral()> SequentialSearch(std::vector< std::function< BooleanOrIntegerLiteral()> > heuristics)
OR-Tools root namespace.
BlossomGraph::NodeIndex NodeIndex
const bool DEBUG_MODE
Definition radix_sort.h:266
std::string FormatCounter(int64_t num)
Definition logging.cc:30
static IntegerLiteral GreaterOrEqual(IntegerVariable i, IntegerValue bound)
static IntegerLiteral LowerOrEqual(IntegerVariable i, IntegerValue bound)
std::function< BooleanOrIntegerLiteral()> fixed_search