Google OR-Tools v9.14
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"
43#include "ortools/sat/sat_parameters.pb.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 = sat_solver_->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 = sat_solver_->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_->RestoreSolverToAssumptionLevel();
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_->RestoreSolverToAssumptionLevel()) {
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) continue;
607 }
608
609 if (sat_solver_->CurrentDecisionLevel() == 0) {
610 if (!LevelZeroLogic()) {
611 return sat_solver_->UnsatStatus();
612 }
613 }
614
615 // Dive: Follow the branch with lowest objective.
616 // Note that we do not creates new nodes here.
617 //
618 // TODO(user): If we have new information and our current objective bound
619 // is higher than any bound in a whole subtree, we might want to just
620 // restart this subtree exploration?
621 while (true) {
622 const int size = current_branch_.size();
623 const int level = sat_solver_->CurrentDecisionLevel();
624
625 // Invariant are tricky:
626 // current_branch_ contains one entry per decision taken + the last one
627 // which we are about to take. If we don't have the last entry, it means
628 // we are about to take a new decision.
629 DCHECK(size == level || size == level + 1);
630 if (size == level) break; // Take new decision.
631
632 const NodeIndex node_index = current_branch_[level];
633 Node& node = nodes_[node_index];
634 DCHECK_GT(node.true_child, node_index);
635 DCHECK_GT(node.false_child, node_index);
636
637 // If the bound of this node is high, restart the main loop..
638 node.UpdateObjective(std::max(
639 current_objective_lb_, integer_trail_->LowerBound(objective_var_)));
640 if (node.MinObjective() > current_objective_lb_) break;
641 DCHECK_EQ(node.MinObjective(), current_objective_lb_) << level;
642
643 // This will be set to the next node index.
644 NodeIndex n;
645 DCHECK(!node.is_deleted);
646 const Literal node_literal = node.Decision();
647
648 // If the variable is already fixed, we bypass the node and connect
649 // its parent directly to the relevant child.
650 if (assignment_.LiteralIsAssigned(node_literal)) {
651 IntegerValue new_lb;
652 if (assignment_.LiteralIsTrue(node_literal)) {
653 n = node.true_child;
654 new_lb = node.true_objective;
655 } else {
656 n = node.false_child;
657 new_lb = node.false_objective;
658 }
659 MarkAsDeletedNodeAndUnreachableSubtree(node);
660
661 // We jump directly to the subnode.
662 // Else we will change the root.
663 current_branch_.pop_back();
664 if (!current_branch_.empty()) {
665 const NodeIndex parent = current_branch_.back();
666 DCHECK(!nodes_[parent].is_deleted);
667 const Literal parent_literal = nodes_[parent].Decision();
668 if (assignment_.LiteralIsTrue(parent_literal)) {
669 nodes_[parent].true_child = n;
670 nodes_[parent].UpdateTrueObjective(new_lb);
671 } else {
672 DCHECK(assignment_.LiteralIsFalse(parent_literal));
673 nodes_[parent].false_child = n;
674 nodes_[parent].UpdateFalseObjective(new_lb);
675 }
676 if (new_lb > current_objective_lb_) {
677 // This is probably not needed.
678 if (n < nodes_.size() && !nodes_[n].IsLeaf()) {
679 current_branch_.push_back(n);
680 nodes_[n].UpdateObjective(new_lb);
681 }
682 break;
683 }
684 } else {
685 if (n >= nodes_.size()) {
686 // We never explored the other branch, so we can just clear all
687 // nodes.
688 num_nodes_in_tree_ = 0;
689 nodes_.clear();
690 } else if (nodes_[n].IsLeaf()) {
691 // Keep the saved basis.
692 num_nodes_in_tree_ = 1;
693 Node root = std::move(nodes_[n]);
694 nodes_.clear();
695 nodes_.push_back(std::move(root));
696 n = NodeIndex(0);
697 } else {
698 // We always make sure the root is at zero.
699 // The root is no longer at zero, that might cause issue.
700 // Cleanup.
701 }
702 }
703 } else {
704 // See if we have better bounds using the current LP state.
705 ExploitReducedCosts(current_branch_[level]);
706 if (node.MinObjective() > current_objective_lb_) break;
707
708 // If both lower bound are the same, we pick the literal branch. We do
709 // that because this is the polarity that was chosen by the SAT
710 // heuristic in the first place. We tried random, it doesn't seems to
711 // work as well.
712 num_decisions_taken_++;
713 const bool choose_true = node.true_objective <= node.false_objective;
714 Literal next_decision;
715 if (choose_true) {
716 n = node.true_child;
717 next_decision = node_literal;
718 } else {
719 n = node.false_child;
720 next_decision = node_literal.Negated();
721 }
722
723 // If we are taking this branch for the first time, we enable the LP and
724 // make sure we solve it before taking the decision. This allow to have
725 // proper pseudo-costs, and also be incremental for the decision we are
726 // about to take.
727 //
728 // We also enable the LP if we have no basis info for this node.
729 if (SaveLpBasisOption() &&
730 (n >= nodes_.size() || !NodeHasBasis(node))) {
731 const auto cleanup =
732 absl::MakeCleanup(UpdateLpIters(&num_lp_iters_save_basis_));
733
734 VLOG(1) << "~~~~";
735 EnableLpAndLoadBestBasis();
736 const int level = sat_solver_->CurrentDecisionLevel();
737 if (!sat_solver_->FinishPropagation()) {
738 return sat_solver_->UnsatStatus();
739 }
740 if (sat_solver_->CurrentDecisionLevel() < level) {
741 node.UpdateObjective(kMaxIntegerValue);
742 break;
743 }
744
745 // The decision might have become assigned, in which case we loop.
746 if (assignment_.LiteralIsAssigned(next_decision)) {
747 continue;
748 }
749
750 SaveLpBasisInto(node);
751
752 // If we are not at the end, disable the LP propagation.
753 if (n < nodes_.size()) {
754 lp_constraint_->EnablePropagation(false);
755 }
756 }
757
758 // Take the decision.
759 const auto cleanup =
760 absl::MakeCleanup(UpdateLpIters(&num_lp_iters_first_branch_));
761 DCHECK(!assignment_.LiteralIsAssigned(next_decision));
762 search_helper_->TakeDecision(next_decision);
763
764 // Conflict?
765 if (current_branch_.size() != sat_solver_->CurrentDecisionLevel()) {
766 MarkBranchAsInfeasible(node, choose_true);
767 break;
768 }
769
770 // Update the proper field and abort the dive if we crossed the
771 // threshold.
772 const IntegerValue lb = integer_trail_->LowerBound(objective_var_);
773 if (choose_true) {
774 node.UpdateTrueObjective(lb);
775 } else {
776 node.UpdateFalseObjective(lb);
777 }
778
779 if (n < nodes_.size()) {
780 nodes_[n].UpdateObjective(lb);
781 } else if (SaveLpBasisOption()) {
782 SaveLpBasisInto(nodes_[CreateNewEmptyNodeIfNeeded()]);
783 }
784
785 if (lb > current_objective_lb_) break;
786 }
787
788 if (VLOG_IS_ON(1)) {
789 shared_response_->LogMessageWithThrottling(
790 "TreeS", absl::StrCat(" (", SmallProgressString(), ")"));
791 }
792
793 if (n < nodes_.size() && !nodes_[n].IsLeaf()) {
794 current_branch_.push_back(n);
795 } else {
796 break;
797 }
798 }
799
800 // If a conflict occurred, we will backtrack.
801 if (current_branch_.size() != sat_solver_->CurrentDecisionLevel()) {
802 continue;
803 }
804
805 // TODO(user): The code is hard to follow. Fix and merge that with test
806 // below.
807 if (!current_branch_.empty()) {
808 const Node& final_node = nodes_[current_branch_.back()];
809 if (assignment_.LiteralIsTrue(final_node.Decision())) {
810 if (final_node.true_objective > current_objective_lb_) {
811 continue;
812 }
813 } else {
814 DCHECK(assignment_.LiteralIsFalse(final_node.Decision()));
815 if (final_node.false_objective > current_objective_lb_) {
816 continue;
817 }
818 }
819 }
820
821 // This test allow to not take a decision when the branch is already closed
822 // (i.e. the true branch or false branch lb is high enough). Adding it
823 // basically changes if we take the decision later when we explore the
824 // branch or right now.
825 //
826 // I feel taking it later is better. It also avoid creating unneeded nodes.
827 // It does change the behavior on a few problem though. For instance on
828 // irp.mps.gz, the search works better without this, whatever the random
829 // seed. Not sure why, maybe it creates more diversity?
830 //
831 // Another difference is that if the search is done and we have a feasible
832 // solution, we will not report it because of this test (except if we are
833 // at the optimal).
834 if (integer_trail_->LowerBound(objective_var_) > current_objective_lb_) {
835 continue;
836 }
837
838 const auto cleanup = absl::MakeCleanup(UpdateLpIters(&num_lp_iters_dive_));
839
840 if (current_branch_.empty()) {
841 VLOG(2) << "DIVE from empty tree";
842 } else {
843 VLOG(2) << "DIVE from " << NodeDebugString(current_branch_.back());
844 }
845
846 if (SaveLpBasisOption() && !lp_constraint_->PropagationIsEnabled()) {
847 // This reuse or create a node to store the basis.
848 const NodeIndex index = CreateNewEmptyNodeIfNeeded();
849
850 EnableLpAndLoadBestBasis();
851 const int level = sat_solver_->CurrentDecisionLevel();
852 if (!sat_solver_->FinishPropagation()) {
853 return sat_solver_->UnsatStatus();
854 }
855
856 // Loop on backtrack or bound improvement.
857 if (sat_solver_->CurrentDecisionLevel() < level) {
858 Node& node = nodes_[index];
859 node.UpdateObjective(kMaxIntegerValue);
860 continue;
861 }
862
863 SaveLpBasisInto(nodes_[index]);
864
865 const IntegerValue obj_lb = integer_trail_->LowerBound(objective_var_);
866 if (obj_lb > current_objective_lb_) {
867 nodes_[index].UpdateObjective(obj_lb);
868 if (!current_branch_.empty()) {
869 Node& parent_node = nodes_[current_branch_.back()];
870 const Literal node_literal = parent_node.Decision();
871 DCHECK(assignment_.LiteralIsAssigned(node_literal));
872 if (assignment_.LiteralIsTrue(node_literal)) {
873 parent_node.UpdateTrueObjective(obj_lb);
874 } else {
875 parent_node.UpdateFalseObjective(obj_lb);
876 }
877 }
878 continue;
879 }
880 }
881
882 // Invariant: The current branch is fully assigned, and the solver is in
883 // sync. And we are not on a "bad" path.
884 const int base_level = sat_solver_->CurrentDecisionLevel();
885 if (DEBUG_MODE) {
886 CHECK_EQ(base_level, current_branch_.size());
887 for (const NodeIndex index : current_branch_) {
888 CHECK(!nodes_[index].is_deleted);
889 const Literal decision = nodes_[index].Decision();
890 if (assignment_.LiteralIsTrue(decision)) {
891 CHECK_EQ(nodes_[index].true_objective, current_objective_lb_);
892 } else {
893 CHECK(assignment_.LiteralIsFalse(decision));
894 CHECK_EQ(nodes_[index].false_objective, current_objective_lb_);
895 }
896 }
897 }
898
899 // We are about to take a new decision, what we will do is dive until
900 // the objective lower bound increase. we will then create a bunch of new
901 // nodes in the tree.
902 //
903 // By analyzing the reason for the increase, we can create less nodes than
904 // if we just followed the initial heuristic.
905 //
906 // TODO(user): In multithread, this change the behavior a lot since we
907 // dive until we beat the best shared bound. Maybe we shouldn't do that.
908 while (true) {
909 // TODO(user): We sometimes branch on the objective variable, this should
910 // probably be avoided.
911 if (sat_solver_->ModelIsUnsat()) return sat_solver_->UnsatStatus();
912 LiteralIndex decision;
913 if (!search_helper_->GetDecision(search_heuristic_, &decision)) {
914 continue;
915 }
916
917 // No new decision: search done.
918 if (time_limit_->LimitReached()) return SatSolver::LIMIT_REACHED;
919 if (decision == kNoLiteralIndex) {
920 feasible_solution_observer();
921 break;
922 }
923
924 num_decisions_taken_++;
925 if (!search_helper_->TakeDecision(Literal(decision))) {
926 return sat_solver_->UnsatStatus();
927 }
928 if (sat_solver_->CurrentDecisionLevel() < base_level) {
929 // TODO(user): it would be nice to mark some node as infeasible if
930 // this is the case. However this could happen after many decision and
931 // we realize with the lp that one of them should have been fixed
932 // earlier, without any infeasibility in the current branch.
933 break;
934 }
935 if (integer_trail_->LowerBound(objective_var_) > current_objective_lb_) {
936 break;
937 }
938 }
939
940 if (sat_solver_->CurrentDecisionLevel() <= base_level) continue;
941
942 // Analyse the reason for objective increase. Deduce a set of new nodes to
943 // append to the tree.
944 //
945 // TODO(user): Try to minimize the number of decisions?
946 const std::vector<Literal> reason =
947 integer_trail_->ReasonFor(IntegerLiteral::GreaterOrEqual(
948 objective_var_, integer_trail_->LowerBound(objective_var_)));
949 std::vector<Literal> decisions = ExtractDecisions(base_level, reason);
950
951 // Bump activities.
952 sat_decision_->BumpVariableActivities(reason);
953 sat_decision_->BumpVariableActivities(decisions);
954 sat_decision_->UpdateVariableActivityIncrement();
955
956 // Create one node per new decisions.
957 DCHECK_EQ(current_branch_.size(), base_level);
958 for (const Literal d : decisions) {
959 AppendNewNodeToCurrentBranch(d);
960 }
961
962 // TODO(user): We should probably save the basis in more cases.
963 if (SaveLpBasisOption() && decisions.size() == 1) {
964 SaveLpBasisInto(nodes_[CreateNewEmptyNodeIfNeeded()]);
965 }
966
967 // Update the objective of the last node in the branch since we just
968 // improved that.
969 if (!current_branch_.empty()) {
970 Node& n = nodes_[current_branch_.back()];
971 if (assignment_.LiteralIsTrue(n.Decision())) {
972 n.UpdateTrueObjective(integer_trail_->LowerBound(objective_var_));
973 } else {
974 n.UpdateFalseObjective(integer_trail_->LowerBound(objective_var_));
975 }
976 }
977
978 // Reset the solver to a correct state since we have a subset of the
979 // current propagation. We backtrack as little as possible.
980 //
981 // The decision level is the number of decision taken.
982 // Decision()[level] is the decision at that level.
983 int backtrack_level = base_level;
984 DCHECK_LE(current_branch_.size(), sat_solver_->CurrentDecisionLevel());
985 while (backtrack_level < current_branch_.size() &&
986 sat_solver_->Decisions()[backtrack_level].literal.Index() ==
987 nodes_[current_branch_[backtrack_level]].literal_index) {
988 ++backtrack_level;
989 }
990 sat_solver_->Backtrack(backtrack_level);
991
992 // Update bounds with reduced costs info.
993 //
994 // TODO(user): Uses old optimal constraint that we just potentially
995 // backtracked over?
996 //
997 // TODO(user): We could do all at once rather than in O(#decision * #size).
998 for (int i = backtrack_level; i < current_branch_.size(); ++i) {
999 ExploitReducedCosts(current_branch_[i]);
1000 }
1001 }
1002
1004}
1005
1006std::vector<Literal> LbTreeSearch::ExtractDecisions(
1007 int base_level, absl::Span<const Literal> conflict) {
1008 std::vector<int> num_per_level(sat_solver_->CurrentDecisionLevel() + 1, 0);
1009 std::vector<bool> is_marked;
1010 for (const Literal l : conflict) {
1011 const AssignmentInfo& info = trail_->Info(l.Variable());
1012 if (info.level <= base_level) continue;
1013 num_per_level[info.level]++;
1014 if (info.trail_index >= is_marked.size()) {
1015 is_marked.resize(info.trail_index + 1);
1016 }
1017 is_marked[info.trail_index] = true;
1018 }
1019
1020 std::vector<Literal> result;
1021 if (is_marked.empty()) return result;
1022 for (int i = is_marked.size() - 1; i >= 0; --i) {
1023 if (!is_marked[i]) continue;
1024
1025 const Literal l = (*trail_)[i];
1026 const AssignmentInfo& info = trail_->Info(l.Variable());
1027 if (info.level <= base_level) break;
1028 if (num_per_level[info.level] == 1) {
1029 result.push_back(l);
1030 continue;
1031 }
1032
1033 // Expand.
1034 num_per_level[info.level]--;
1035 for (const Literal new_l : trail_->Reason(l.Variable())) {
1036 const AssignmentInfo& new_info = trail_->Info(new_l.Variable());
1037 if (new_info.level <= base_level) continue;
1038 if (is_marked[new_info.trail_index]) continue;
1039 is_marked[new_info.trail_index] = true;
1040 num_per_level[new_info.level]++;
1041 }
1042 }
1043
1044 // We prefer to keep the same order.
1045 std::reverse(result.begin(), result.end());
1046 return result;
1047}
1048
1049LbTreeSearch::NodeIndex LbTreeSearch::CreateNewEmptyNodeIfNeeded() {
1050 NodeIndex n(0);
1051 if (current_branch_.empty()) {
1052 if (nodes_.empty()) {
1053 ++num_nodes_in_tree_;
1054 nodes_.emplace_back(current_objective_lb_);
1055 } else {
1056 DCHECK_EQ(nodes_.size(), 1);
1057 }
1058 } else {
1059 const NodeIndex parent = current_branch_.back();
1060 DCHECK(!nodes_[parent].is_deleted);
1061 const Literal parent_literal = nodes_[parent].Decision();
1062 if (assignment_.LiteralIsTrue(parent_literal)) {
1063 if (nodes_[parent].true_child >= nodes_.size()) {
1064 n = NodeIndex(nodes_.size());
1065 ++num_nodes_in_tree_;
1066 nodes_[parent].true_child = NodeIndex(nodes_.size());
1067 nodes_.emplace_back(current_objective_lb_);
1068 } else {
1069 n = nodes_[parent].true_child;
1070 }
1071 nodes_[parent].UpdateTrueObjective(current_objective_lb_);
1072 } else {
1073 DCHECK(assignment_.LiteralIsFalse(parent_literal));
1074 if (nodes_[parent].false_child >= nodes_.size()) {
1075 n = NodeIndex(nodes_.size());
1076 ++num_nodes_in_tree_;
1077 nodes_[parent].false_child = NodeIndex(nodes_.size());
1078 nodes_.emplace_back(current_objective_lb_);
1079 } else {
1080 n = nodes_[parent].false_child;
1081 }
1082 nodes_[parent].UpdateFalseObjective(current_objective_lb_);
1083 }
1084 }
1085 DCHECK(!nodes_[n].is_deleted);
1086 DCHECK_EQ(nodes_[n].literal_index, kNoLiteralIndex);
1087 return n;
1088}
1089
1090void LbTreeSearch::AppendNewNodeToCurrentBranch(Literal decision) {
1091 NodeIndex n(nodes_.size());
1092 if (current_branch_.empty()) {
1093 if (nodes_.empty()) {
1094 ++num_nodes_in_tree_;
1095 nodes_.emplace_back(current_objective_lb_);
1096 } else {
1097 DCHECK_EQ(nodes_.size(), 1);
1098 n = 0;
1099 }
1100 } else {
1101 const NodeIndex parent = current_branch_.back();
1102 DCHECK(!nodes_[parent].is_deleted);
1103 const Literal parent_literal = nodes_[parent].Decision();
1104 if (assignment_.LiteralIsTrue(parent_literal)) {
1105 if (nodes_[parent].true_child < nodes_.size()) {
1106 n = nodes_[parent].true_child;
1107 } else {
1108 ++num_nodes_in_tree_;
1109 nodes_.emplace_back(current_objective_lb_);
1110 nodes_[parent].true_child = n;
1111 }
1112 nodes_[parent].UpdateTrueObjective(current_objective_lb_);
1113 } else {
1114 DCHECK(assignment_.LiteralIsFalse(parent_literal));
1115 if (nodes_[parent].false_child < nodes_.size()) {
1116 n = nodes_[parent].false_child;
1117 } else {
1118 ++num_nodes_in_tree_;
1119 nodes_.emplace_back(current_objective_lb_);
1120 nodes_[parent].false_child = n;
1121 }
1122 nodes_[parent].UpdateFalseObjective(current_objective_lb_);
1123 }
1124 }
1125 DCHECK_LT(n, nodes_.size());
1126 DCHECK_EQ(nodes_[n].literal_index, kNoLiteralIndex) << " issue " << n;
1127 nodes_[n].SetDecision(decision);
1128 nodes_[n].UpdateObjective(current_objective_lb_);
1129 current_branch_.push_back(n);
1130}
1131
1132// Looking at the reduced costs, we can already have a bound for one of the
1133// branch. Increasing the corresponding objective can save some branches,
1134// and also allow for a more incremental LP solving since we do less back
1135// and forth.
1136//
1137// TODO(user): The code to recover that is a bit convoluted. Alternatively
1138// Maybe we should do a "fast" propagation without the LP in each branch.
1139// That will work as long as we keep these optimal LP constraints around
1140// and propagate them.
1141//
1142// TODO(user): Incorporate this in the heuristic so we choose more Booleans
1143// inside these LP explanations?
1144void LbTreeSearch::ExploitReducedCosts(NodeIndex n) {
1145 if (lp_constraint_ == nullptr) return;
1146
1147 // TODO(user): we could consider earlier constraints instead of just
1148 // looking at the last one, but experiments didn't really show a big
1149 // gain.
1150 const auto& cts = lp_constraint_->OptimalConstraints();
1151 if (cts.empty()) return;
1152 const std::unique_ptr<IntegerSumLE128>& rc = cts.back();
1153
1154 // Note that this return literal EQUIVALENT to the node.literal, not just
1155 // implied by it. We need that for correctness.
1156 int num_tests = 0;
1157 Node& node = nodes_[n];
1158 DCHECK(!node.is_deleted);
1159 const Literal node_literal = node.Decision();
1160 DCHECK(!assignment_.LiteralIsAssigned(node_literal));
1161 for (const IntegerLiteral integer_literal :
1162 integer_encoder_->GetIntegerLiterals(node_literal)) {
1163 // To avoid bad corner case. Not sure it ever triggers.
1164 if (++num_tests > 10) break;
1165
1166 const std::pair<IntegerValue, IntegerValue> bounds =
1167 rc->ConditionalLb(integer_literal, objective_var_);
1168 if (bounds.first > node.false_objective) {
1169 ++num_rc_detected_;
1170 node.UpdateFalseObjective(bounds.first);
1171 }
1172 if (bounds.second > node.true_objective) {
1173 ++num_rc_detected_;
1174 node.UpdateTrueObjective(bounds.second);
1175 }
1176 }
1177}
1178
1179} // namespace sat
1180} // namespace operations_research
An helper class to share the code used by the different kind of search.
SatSolver::Status Search(const std::function< void()> &feasible_solution_observer)
Explores the search space.
A class that stores the collection of all LP constraints in a model.
int64_t num_lp_changes() const
This can serve as a timestamp to know if a saved basis is out of date.
T Get(std::function< T(const Model &)> f) const
Similar to Add() but this is const.
Definition model.h:93
const AssignmentInfo & Info(BooleanVariable var) const
Definition sat_base.h:463
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::string FormatCounter(int64_t num)
Prints a positive number with separators for easier reading (ex: 1'348'065).
Definition util.cc:44
std::function< BooleanOrIntegerLiteral()> SatSolverHeuristic(Model *model)
Returns the BooleanOrIntegerLiteral advised by the underlying SAT solver.
std::function< BooleanOrIntegerLiteral()> LpPseudoCostHeuristic(Model *model)
std::function< BooleanOrIntegerLiteral()> MostFractionalHeuristic(Model *model)
Choose the variable with most fractional LP value.
std::function< BooleanOrIntegerLiteral()> SequentialSearch(std::vector< std::function< BooleanOrIntegerLiteral()> > heuristics)
In SWIG mode, we don't want anything besides these top-level includes.
BlossomGraph::NodeIndex NodeIndex
const bool DEBUG_MODE
Definition radix_sort.h:266
Information about a variable assignment.
Definition sat_base.h:242
int32_t trail_index
The index of this assignment in the trail.
Definition sat_base.h:263
static IntegerLiteral GreaterOrEqual(IntegerVariable i, IntegerValue bound)
static IntegerLiteral LowerOrEqual(IntegerVariable i, IntegerValue bound)
std::function< BooleanOrIntegerLiteral()> fixed_search
Fixed search, built from above building blocks.