OR-Tools  8.2
perfect_matching.cc
Go to the documentation of this file.
1// Copyright 2010-2018 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 "absl/memory/memory.h"
18
19namespace operations_research {
20
21void MinCostPerfectMatching::Reset(int num_nodes) {
22 graph_ = absl::make_unique<BlossomGraph>(num_nodes);
23 optimal_cost_ = 0;
24 matches_.assign(num_nodes, -1);
25}
26
28 CHECK_GE(cost, 0) << "Not supported for now, just shift your costs.";
29 if (tail == head) {
30 VLOG(1) << "Ignoring self-arc: " << tail << " <-> " << head
31 << " cost: " << cost;
32 return;
33 }
34 maximum_edge_cost_ = std::max(maximum_edge_cost_, cost);
37}
38
40 optimal_solution_found_ = false;
41
42 // We want all dual and all slack value to never overflow. After Initialize()
43 // they are both bounded by the 2 * maximum cost. And we track an upper bound
44 // on these quantities. The factor two is because of the re-scaling we do
45 // internally since all our dual values are actually multiple of 1/2.
46 //
47 // Note that since the whole code in BlossomGraph assumes that dual/slack have
48 // a magnitude that is always lower than kMaxCostValue it is important to use
49 // it here since there is no reason it cannot be smaller than kint64max.
50 //
51 // TODO(user): Improve the overflow detection if needed. The current one seems
52 // ok though.
53 int64 overflow_detection = CapAdd(maximum_edge_cost_, maximum_edge_cost_);
54 if (overflow_detection >= BlossomGraph::kMaxCostValue) {
55 return Status::INTEGER_OVERFLOW;
56 }
57
58 const int num_nodes = matches_.size();
59 if (!graph_->Initialize()) return Status::INFEASIBLE;
60 VLOG(2) << graph_->DebugString();
61 VLOG(1) << "num_unmatched: " << num_nodes - graph_->NumMatched()
62 << " dual_objective: " << graph_->DualObjective();
63
64 while (graph_->NumMatched() != num_nodes) {
65 graph_->PrimalUpdates();
66 if (DEBUG_MODE) {
67 graph_->DebugCheckNoPossiblePrimalUpdates();
68 }
69
70 VLOG(1) << "num_unmatched: " << num_nodes - graph_->NumMatched()
71 << " dual_objective: " << graph_->DualObjective();
72 if (graph_->NumMatched() == num_nodes) break;
73
75 graph_->ComputeMaxCommonTreeDualDeltaAndResetPrimalEdgeQueue();
76 overflow_detection = CapAdd(overflow_detection, std::abs(delta.value()));
77 if (overflow_detection >= BlossomGraph::kMaxCostValue) {
78 return Status::INTEGER_OVERFLOW;
79 }
80
81 if (delta == 0) break; // Infeasible!
82 graph_->UpdateAllTrees(delta);
83 }
84
85 VLOG(1) << "End: " << graph_->NumMatched() << " / " << num_nodes;
86 graph_->DisplayStats();
87 if (graph_->NumMatched() < num_nodes) {
88 return Status::INFEASIBLE;
89 }
90 VLOG(2) << graph_->DebugString();
91 CHECK(graph_->DebugDualsAreFeasible());
92
93 // TODO(user): Maybe there is a faster/better way to recover the mapping
94 // in the presence of blossoms.
95 graph_->ExpandAllBlossoms();
96 for (int i = 0; i < num_nodes; ++i) {
97 matches_[i] = graph_->Match(BlossomGraph::NodeIndex(i)).value();
98 }
99
100 optimal_solution_found_ = true;
101 optimal_cost_ = graph_->DualObjective().value();
102 if (optimal_cost_ == kint64max) return Status::COST_OVERFLOW;
103 return Status::OPTIMAL;
104}
105
108
111const BlossomGraph::EdgeIndex BlossomGraph::kNoEdgeIndex =
112 BlossomGraph::EdgeIndex(-1);
115
117 graph_.resize(num_nodes);
118 nodes_.reserve(num_nodes);
119 root_blossom_node_.resize(num_nodes);
120 for (NodeIndex n(0); n < num_nodes; ++n) {
121 root_blossom_node_[n] = n;
122 nodes_.push_back(Node(n));
123 }
124}
125
127 DCHECK_GE(tail, 0);
128 DCHECK_LT(tail, nodes_.size());
129 DCHECK_GE(head, 0);
130 DCHECK_LT(head, nodes_.size());
131 DCHECK_GE(cost, 0);
132 DCHECK(!is_initialized_);
133 const EdgeIndex index(edges_.size());
134 edges_.push_back(Edge(tail, head, cost));
135 graph_[tail].push_back(index);
136 graph_[head].push_back(index);
137}
138
139// TODO(user): Code the more advanced "Fractional matching initialization"
140// heuristic.
141//
142// TODO(user): Add a preprocessing step that performs the 'forced' matches?
144 CHECK(!is_initialized_);
145 is_initialized_ = true;
146
147 for (NodeIndex n(0); n < nodes_.size(); ++n) {
148 if (graph_[n].empty()) return false; // INFEASIBLE.
149 CostValue min_cost = kMaxCostValue;
150
151 // Initialize the dual of each nodes to min_cost / 2.
152 //
153 // TODO(user): We might be able to do better for odd min_cost, but then
154 // we might need to scale by 4? think about it.
155 for (const EdgeIndex e : graph_[n]) {
156 min_cost = std::min(min_cost, edges_[e].pseudo_slack);
157 }
158 DCHECK_NE(min_cost, kMaxCostValue);
159 nodes_[n].pseudo_dual = min_cost / 2;
160
161 // Starts with all nodes as tree roots.
162 nodes_[n].type = 1;
163 }
164
165 // Update the slack of each edges now that nodes might have non-zero duals.
166 // Note that we made sure that all updated slacks are non-negative.
167 for (EdgeIndex e(0); e < edges_.size(); ++e) {
168 Edge& mutable_edge = edges_[e];
169 mutable_edge.pseudo_slack -= nodes_[mutable_edge.tail].pseudo_dual +
170 nodes_[mutable_edge.head].pseudo_dual;
171 DCHECK_GE(mutable_edge.pseudo_slack, 0);
172 }
173
174 for (NodeIndex n(0); n < nodes_.size(); ++n) {
175 if (NodeIsMatched(n)) continue;
176
177 // After this greedy update, there will be at least an edge with a
178 // slack of zero.
179 CostValue min_slack = kMaxCostValue;
180 for (const EdgeIndex e : graph_[n]) {
181 min_slack = std::min(min_slack, edges_[e].pseudo_slack);
182 }
183 DCHECK_NE(min_slack, kMaxCostValue);
184 if (min_slack > 0) {
185 nodes_[n].pseudo_dual += min_slack;
186 for (const EdgeIndex e : graph_[n]) {
187 edges_[e].pseudo_slack -= min_slack;
188 }
189 DebugUpdateNodeDual(n, min_slack);
190 }
191
192 // Match this node if possible.
193 //
194 // TODO(user): Optimize by merging this loop with the one above?
195 for (const EdgeIndex e : graph_[n]) {
196 const Edge& edge = edges_[e];
197 if (edge.pseudo_slack != 0) continue;
198 if (!NodeIsMatched(edge.OtherEnd(n))) {
199 nodes_[edge.tail].type = 0;
200 nodes_[edge.tail].match = edge.head;
201 nodes_[edge.head].type = 0;
202 nodes_[edge.head].match = edge.tail;
203 break;
204 }
205 }
206 }
207
208 // Initialize unmatched_nodes_.
209 for (NodeIndex n(0); n < nodes_.size(); ++n) {
210 if (NodeIsMatched(n)) continue;
211 unmatched_nodes_.push_back(n);
212 }
213
214 // Scale everything by 2 and update the dual cost. Note that we made sure that
215 // there cannot be an integer overflow at the beginning of Solve().
216 //
217 // This scaling allows to only have integer weights during the algorithm
218 // because the slack of [+] -- [+] edges will always stay even.
219 //
220 // TODO(user): Reduce the number of loops we do in the initialization. We
221 // could likely just scale the edge cost as we fill them.
222 for (NodeIndex n(0); n < nodes_.size(); ++n) {
223 DCHECK_LE(nodes_[n].pseudo_dual, kMaxCostValue / 2);
224 nodes_[n].pseudo_dual *= 2;
225 AddToDualObjective(nodes_[n].pseudo_dual);
226#ifndef NDEBUG
227 nodes_[n].dual = nodes_[n].pseudo_dual;
228#endif
229 }
230 for (EdgeIndex e(0); e < edges_.size(); ++e) {
231 DCHECK_LE(edges_[e].pseudo_slack, kMaxCostValue / 2);
232 edges_[e].pseudo_slack *= 2;
233#ifndef NDEBUG
234 edges_[e].slack = edges_[e].pseudo_slack;
235#endif
236 }
237
238 // Initialize the edge priority queues and the primal update queue.
239 // We only need to do that if we have unmatched nodes.
240 if (!unmatched_nodes_.empty()) {
241 primal_update_edge_queue_.clear();
242 for (EdgeIndex e(0); e < edges_.size(); ++e) {
243 Edge& edge = edges_[e];
244 const bool tail_is_plus = nodes_[edge.tail].IsPlus();
245 const bool head_is_plus = nodes_[edge.head].IsPlus();
246 if (tail_is_plus && head_is_plus) {
247 plus_plus_pq_.Add(&edge);
248 if (edge.pseudo_slack == 0) primal_update_edge_queue_.push_back(e);
249 } else if (tail_is_plus || head_is_plus) {
250 plus_free_pq_.Add(&edge);
251 if (edge.pseudo_slack == 0) primal_update_edge_queue_.push_back(e);
252 }
253 }
254 }
255
256 return true;
257}
258
260 // TODO(user): Avoid this linear loop.
261 CostValue best_update = kMaxCostValue;
262 for (NodeIndex n(0); n < nodes_.size(); ++n) {
263 const Node& node = nodes_[n];
264 if (node.IsBlossom() && node.IsMinus()) {
265 best_update = std::min(best_update, Dual(node));
266 }
267 }
268
269 // This code only works because all tree_dual_delta are the same.
270 CHECK(!unmatched_nodes_.empty());
271 const CostValue tree_delta = nodes_[unmatched_nodes_.front()].tree_dual_delta;
272 CostValue plus_plus_slack = kMaxCostValue;
273 if (!plus_plus_pq_.IsEmpty()) {
274 DCHECK_EQ(plus_plus_pq_.Top()->pseudo_slack % 2, 0) << "Non integer bound!";
275 plus_plus_slack = plus_plus_pq_.Top()->pseudo_slack / 2 - tree_delta;
276 best_update = std::min(best_update, plus_plus_slack);
277 }
278 CostValue plus_free_slack = kMaxCostValue;
279 if (!plus_free_pq_.IsEmpty()) {
280 plus_free_slack = plus_free_pq_.Top()->pseudo_slack - tree_delta;
281 best_update = std::min(best_update, plus_free_slack);
282 }
283
284 // This means infeasible, and returning zero will abort the search.
285 if (best_update == kMaxCostValue) return CostValue(0);
286
287 // Initialize primal_update_edge_queue_ with all the edges that will have a
288 // slack of zero once we apply the update.
289 //
290 // NOTE(user): If we want more "determinism" and be independent on the PQ
291 // algorithm, we could std::sort() the primal_update_edge_queue_ here.
292 primal_update_edge_queue_.clear();
293 if (plus_plus_slack == best_update) {
294 plus_plus_pq_.AllTop(&tmp_all_tops_);
295 for (const Edge* pt : tmp_all_tops_) {
296 primal_update_edge_queue_.push_back(EdgeIndex(pt - &edges_.front()));
297 }
298 }
299 if (plus_free_slack == best_update) {
300 plus_free_pq_.AllTop(&tmp_all_tops_);
301 for (const Edge* pt : tmp_all_tops_) {
302 primal_update_edge_queue_.push_back(EdgeIndex(pt - &edges_.front()));
303 }
304 }
305
306 return best_update;
307}
308
310 ++num_dual_updates_;
311
312 // Reminder: the tree roots are exactly the unmatched nodes.
313 CHECK_GE(delta, 0);
314 for (const NodeIndex n : unmatched_nodes_) {
315 CHECK(!NodeIsMatched(n));
316 AddToDualObjective(delta);
317 nodes_[n].tree_dual_delta += delta;
318 }
319
320 if (DEBUG_MODE) {
321 for (NodeIndex n(0); n < nodes_.size(); ++n) {
322 const Node& node = nodes_[n];
323 if (node.IsPlus()) DebugUpdateNodeDual(n, delta);
324 if (node.IsMinus()) DebugUpdateNodeDual(n, -delta);
325 }
326 }
327}
328
330 // An unmatched node must be a tree root.
331 const Node& node = nodes_[n];
332 CHECK(node.match != n || (node.root == n && node.IsPlus()));
333 return node.match != n;
334}
335
337 const Node& node = nodes_[n];
338 if (DEBUG_MODE) {
339 if (node.IsMinus()) CHECK_EQ(node.parent, node.match);
340 if (node.IsPlus()) CHECK_EQ(n, node.match);
341 }
342 return node.match;
343}
344
345// Meant to only be used in DEBUG to make sure our queue in PrimalUpdates()
346// do not miss any potential edges.
348 for (EdgeIndex e(0); e < edges_.size(); ++e) {
349 const Edge& edge = edges_[e];
350 if (Head(edge) == Tail(edge)) continue;
351
352 CHECK(!nodes_[Tail(edge)].is_internal);
353 CHECK(!nodes_[Head(edge)].is_internal);
354 if (Slack(edge) != 0) continue;
355
356 // Make sure tail is a plus node if possible.
357 NodeIndex tail = Tail(edge);
358 NodeIndex head = Head(edge);
359 if (!nodes_[tail].IsPlus()) std::swap(tail, head);
360 if (!nodes_[tail].IsPlus()) continue;
361
362 if (nodes_[head].IsFree()) {
363 VLOG(2) << DebugString();
364 LOG(FATAL) << "Possible Grow! " << tail << " " << head;
365 }
366 if (nodes_[head].IsPlus()) {
367 if (nodes_[tail].root == nodes_[head].root) {
368 LOG(FATAL) << "Possible Shrink!";
369 } else {
370 LOG(FATAL) << "Possible augment!";
371 }
372 }
373 }
374 for (const Node& node : nodes_) {
375 if (node.IsMinus() && node.IsBlossom() && Dual(node) == 0) {
376 LOG(FATAL) << "Possible expand!";
377 }
378 }
379}
380
382 // Any Grow/Augment/Shrink/Expand operation can add new tight edges that need
383 // to be explored again.
384 //
385 // TODO(user): avoid adding duplicates?
386 while (true) {
387 possible_shrink_.clear();
388
389 // First, we Grow/Augment as much as possible.
390 while (!primal_update_edge_queue_.empty()) {
391 const EdgeIndex e = primal_update_edge_queue_.back();
392 primal_update_edge_queue_.pop_back();
393
394 // Because of the Expand() operation, the edge may have become un-tight
395 // since it has been inserted in the tight edges queue. It's cheaper to
396 // detect it here and skip it than it would be to dynamically update the
397 // queue to only keep actually tight edges at all times.
398 const Edge& edge = edges_[e];
399 if (Slack(edge) != 0) continue;
400
401 NodeIndex tail = Tail(edge);
402 NodeIndex head = Head(edge);
403 if (!nodes_[tail].IsPlus()) std::swap(tail, head);
404 if (!nodes_[tail].IsPlus()) continue;
405
406 if (nodes_[head].IsFree()) {
407 Grow(e, tail, head);
408 } else if (nodes_[head].IsPlus()) {
409 if (nodes_[tail].root != nodes_[head].root) {
410 Augment(e);
411 } else {
412 possible_shrink_.push_back(e);
413 }
414 }
415 }
416
417 // Shrink all potential Blossom.
418 for (const EdgeIndex e : possible_shrink_) {
419 const Edge& edge = edges_[e];
420 const NodeIndex tail = Tail(edge);
421 const NodeIndex head = Head(edge);
422 const Node& tail_node = nodes_[tail];
423 const Node& head_node = nodes_[head];
424 if (tail_node.IsPlus() && head_node.IsPlus() &&
425 tail_node.root == head_node.root && tail != head) {
426 Shrink(e);
427 }
428 }
429
430 // Delay expand if any blossom was created.
431 if (!primal_update_edge_queue_.empty()) continue;
432
433 // Expand Blossom if any.
434 //
435 // TODO(user): Avoid doing a O(num_nodes). Also expand all blossom
436 // recursively? I am not sure it is a good heuristic to expand all possible
437 // blossom before trying the other operations though.
438 int num_expands = 0;
439 for (NodeIndex n(0); n < nodes_.size(); ++n) {
440 const Node& node = nodes_[n];
441 if (node.IsMinus() && node.IsBlossom() && Dual(node) == 0) {
442 ++num_expands;
443 Expand(n);
444 }
445 }
446 if (num_expands == 0) break;
447 }
448}
449
451 // The slack of all edge must be non-negative.
452 for (const Edge& edge : edges_) {
453 if (Slack(edge) < 0) return false;
454 }
455
456 // The dual of all Blossom must be non-negative.
457 for (const Node& node : nodes_) {
458 if (node.IsBlossom() && Dual(node) < 0) return false;
459 }
460 return true;
461}
462
464 if (Tail(edge) == Head(edge)) return false;
465 if (nodes_[Tail(edge)].IsInternal()) return false;
466 if (nodes_[Head(edge)].IsInternal()) return false;
467 return Slack(edge) == 0;
468}
469
471 ++num_grows_;
472 VLOG(2) << "Grow " << tail << " -> " << head << " === " << Match(head);
473
475 DCHECK(nodes_[tail].IsPlus());
476 DCHECK(nodes_[head].IsFree());
478
479 const NodeIndex root = nodes_[tail].root;
480 const NodeIndex leaf = Match(head);
481
482 Node& head_node = nodes_[head];
483 head_node.root = root;
484 head_node.parent = tail;
485 head_node.type = -1;
486
487 // head was free and is now a [-] node.
488 const CostValue tree_dual = nodes_[root].tree_dual_delta;
489 head_node.pseudo_dual += tree_dual;
490 for (const NodeIndex subnode : SubNodes(head)) {
491 for (const EdgeIndex e : graph_[subnode]) {
492 Edge& edge = edges_[e];
493 const NodeIndex other_end = OtherEnd(edge, subnode);
494 if (other_end == head) continue;
495 edge.pseudo_slack -= tree_dual;
496 if (plus_free_pq_.Contains(&edge)) plus_free_pq_.Remove(&edge);
497 }
498 }
499
500 Node& leaf_node = nodes_[leaf];
501 leaf_node.root = root;
502 leaf_node.parent = head;
503 leaf_node.type = +1;
504
505 // leaf was free and is now a [+] node.
506 leaf_node.pseudo_dual -= tree_dual;
507 for (const NodeIndex subnode : SubNodes(leaf)) {
508 for (const EdgeIndex e : graph_[subnode]) {
509 Edge& edge = edges_[e];
510 const NodeIndex other_end = OtherEnd(edge, subnode);
511 if (other_end == leaf) continue;
512 edge.pseudo_slack += tree_dual;
513 const Node& other_node = nodes_[other_end];
514 if (other_node.IsPlus()) {
515 // The edge switch from [+] -- [0] to [+] -- [+].
516 DCHECK(plus_free_pq_.Contains(&edge));
517 DCHECK(!plus_plus_pq_.Contains(&edge));
518 plus_free_pq_.Remove(&edge);
519 plus_plus_pq_.Add(&edge);
520 if (edge.pseudo_slack == 2 * tree_dual) {
521 DCHECK_EQ(Slack(edge), 0);
522 primal_update_edge_queue_.push_back(e);
523 }
524 } else if (other_node.IsFree()) {
525 // We have a new [+] -- [0] edge.
526 DCHECK(!plus_free_pq_.Contains(&edge));
527 DCHECK(!plus_plus_pq_.Contains(&edge));
528 plus_free_pq_.Add(&edge);
529 if (edge.pseudo_slack == tree_dual) {
530 DCHECK_EQ(Slack(edge), 0);
531 primal_update_edge_queue_.push_back(e);
532 }
533 }
534 }
535 }
536}
537
538void BlossomGraph::AppendNodePathToRoot(NodeIndex n,
539 std::vector<NodeIndex>* path) const {
540 while (true) {
541 path->push_back(n);
542 n = nodes_[n].parent;
543 if (n == path->back()) break;
544 }
545}
546
547void BlossomGraph::Augment(EdgeIndex e) {
548 ++num_augments_;
549
550 const Edge& edge = edges_[e];
551 VLOG(2) << "Augment " << Tail(edge) << " -> " << Head(edge);
553 DCHECK(nodes_[Tail(edge)].IsPlus());
554 DCHECK(nodes_[Head(edge)].IsPlus());
555
556 const NodeIndex root_a = nodes_[Tail(edge)].root;
557 const NodeIndex root_b = nodes_[Head(edge)].root;
558 DCHECK_NE(root_a, root_b);
559
560 // Compute the path from root_a to root_b.
561 std::vector<NodeIndex> node_path;
562 AppendNodePathToRoot(Tail(edge), &node_path);
563 std::reverse(node_path.begin(), node_path.end());
564 AppendNodePathToRoot(Head(edge), &node_path);
565
566 // TODO(user): Check all dual/slack same after primal op?
567 const CostValue delta_a = nodes_[root_a].tree_dual_delta;
568 const CostValue delta_b = nodes_[root_b].tree_dual_delta;
569 nodes_[root_a].tree_dual_delta = 0;
570 nodes_[root_b].tree_dual_delta = 0;
571
572 // Make all the nodes from both trees free while keeping the
573 // current matching.
574 //
575 // TODO(user): It seems that we may waste some computation since the part of
576 // the tree not in the path between roots can lead to the same Grow()
577 // operations later when one of its node is ratched to a new root.
578 //
579 // TODO(user): Reduce this O(num_nodes) complexity. We might be able to
580 // even do O(num_node_in_path) with lazy updates. Note that this operation
581 // will only be performed at most num_initial_unmatched_nodes / 2 times
582 // though.
583 for (NodeIndex n(0); n < nodes_.size(); ++n) {
584 Node& node = nodes_[n];
585 if (node.IsInternal()) continue;
586 const NodeIndex root = node.root;
587 if (root != root_a && root != root_b) continue;
588
589 const CostValue delta = node.type * (root == root_a ? delta_a : delta_b);
590 node.pseudo_dual += delta;
591 for (const NodeIndex subnode : SubNodes(n)) {
592 for (const EdgeIndex e : graph_[subnode]) {
593 Edge& edge = edges_[e];
594 const NodeIndex other_end = OtherEnd(edge, subnode);
595 if (other_end == n) continue;
596 edge.pseudo_slack -= delta;
597
598 // If the other end is not in one of the two trees, and it is a plus
599 // node, we add it the plus_free queue. All previous [+]--[0] and
600 // [+]--[+] edges need to be removed from the queues.
601 const Node& other_node = nodes_[other_end];
602 if (other_node.root != root_a && other_node.root != root_b &&
603 other_node.IsPlus()) {
604 if (plus_plus_pq_.Contains(&edge)) plus_plus_pq_.Remove(&edge);
605 DCHECK(!plus_free_pq_.Contains(&edge));
606 plus_free_pq_.Add(&edge);
607 if (Slack(edge) == 0) primal_update_edge_queue_.push_back(e);
608 } else {
609 if (plus_plus_pq_.Contains(&edge)) plus_plus_pq_.Remove(&edge);
610 if (plus_free_pq_.Contains(&edge)) plus_free_pq_.Remove(&edge);
611 }
612 }
613 }
614
615 node.type = 0;
616 node.parent = node.root = n;
617 }
618
619 // Change the matching of nodes along node_path.
620 CHECK_EQ(node_path.size() % 2, 0);
621 for (int i = 0; i < node_path.size(); i += 2) {
622 nodes_[node_path[i]].match = node_path[i + 1];
623 nodes_[node_path[i + 1]].match = node_path[i];
624 }
625
626 // Update unmatched_nodes_.
627 //
628 // TODO(user): This could probably be optimized if needed. But we do usually
629 // iterate a lot more over it than we update it. Note that as long as we use
630 // the same delta for all trees, this is not even needed.
631 int new_size = 0;
632 for (const NodeIndex n : unmatched_nodes_) {
633 if (!NodeIsMatched(n)) unmatched_nodes_[new_size++] = n;
634 }
635 CHECK_EQ(unmatched_nodes_.size(), new_size + 2);
636 unmatched_nodes_.resize(new_size);
637}
638
639int BlossomGraph::GetDepth(NodeIndex n) const {
640 int depth = 0;
641 while (true) {
642 const NodeIndex parent = nodes_[n].parent;
643 if (parent == n) break;
644 ++depth;
645 n = parent;
646 }
647 return depth;
648}
649
650void BlossomGraph::Shrink(EdgeIndex e) {
651 ++num_shrinks_;
652
653 const Edge& edge = edges_[e];
655 DCHECK(nodes_[Tail(edge)].IsPlus());
656 DCHECK(nodes_[Head(edge)].IsPlus());
657 DCHECK_EQ(nodes_[Tail(edge)].root, nodes_[Head(edge)].root);
658
659 CHECK_NE(Tail(edge), Head(edge)) << e;
660
661 // Find lowest common ancestor and the two node paths to reach it. Note that
662 // we do not add it to the paths.
663 NodeIndex lca_index = kNoNodeIndex;
664 std::vector<NodeIndex> tail_path;
665 std::vector<NodeIndex> head_path;
666 {
667 NodeIndex tail = Tail(edge);
668 NodeIndex head = Head(edge);
669 int tail_depth = GetDepth(tail);
670 int head_depth = GetDepth(head);
671 if (tail_depth > head_depth) {
672 std::swap(tail, head);
673 std::swap(tail_depth, head_depth);
674 }
675 VLOG(2) << "Shrink " << tail << " <-> " << head;
676
677 while (head_depth > tail_depth) {
678 head_path.push_back(head);
679 head = nodes_[head].parent;
680 --head_depth;
681 }
682 while (tail != head) {
683 DCHECK_EQ(tail_depth, head_depth);
684 DCHECK_GE(tail_depth, 0);
685 if (DEBUG_MODE) {
686 --tail_depth;
687 --head_depth;
688 }
689
690 tail_path.push_back(tail);
691 tail = nodes_[tail].parent;
692
693 head_path.push_back(head);
694 head = nodes_[head].parent;
695 }
696 lca_index = tail;
697 VLOG(2) << "LCA " << lca_index;
698 }
699 Node& lca = nodes_[lca_index];
700 DCHECK(lca.IsPlus());
701
702 // Fill the cycle.
703 std::vector<NodeIndex> blossom = {lca_index};
704 std::reverse(head_path.begin(), head_path.end());
705 blossom.insert(blossom.end(), head_path.begin(), head_path.end());
706 blossom.insert(blossom.end(), tail_path.begin(), tail_path.end());
707 CHECK_EQ(blossom.size() % 2, 1);
708
709 const CostValue tree_dual = nodes_[lca.root].tree_dual_delta;
710
711 // Save all values that will be needed if we expand this Blossom later.
712 CHECK_GT(blossom.size(), 1);
713 Node& backup_node = nodes_[blossom[1]];
714#ifndef NDEBUG
715 backup_node.saved_dual = lca.dual;
716#endif
717 backup_node.saved_pseudo_dual = lca.pseudo_dual + tree_dual;
718
719 // Set the new dual of the node to zero.
720#ifndef NDEBUG
721 lca.dual = 0;
722#endif
723 lca.pseudo_dual = -tree_dual;
724 CHECK_EQ(Dual(lca), 0);
725
726 // Mark node as internal, but do not change their type to zero yet.
727 // We need to do that first to properly detect edges between two internal
728 // nodes in the second loop below.
729 for (const NodeIndex n : blossom) {
730 VLOG(2) << "blossom-node: " << NodeDebugString(n);
731 if (n != lca_index) {
732 nodes_[n].is_internal = true;
733 }
734 }
735
736 // Update the dual of all edges and the priority queueus.
737 for (const NodeIndex n : blossom) {
738 Node& mutable_node = nodes_[n];
739 const bool was_minus = mutable_node.IsMinus();
740 const CostValue slack_adjust =
741 mutable_node.IsMinus() ? tree_dual : -tree_dual;
742 if (n != lca_index) {
743 mutable_node.pseudo_dual -= slack_adjust;
744#ifndef NDEBUG
745 DCHECK_EQ(mutable_node.dual, mutable_node.pseudo_dual);
746#endif
747 mutable_node.type = 0;
748 }
749 for (const NodeIndex subnode : SubNodes(n)) {
750 // Subtle: We update root_blossom_node_ while we loop, so for new internal
751 // edges, depending if an edge "other end" appear after or before, it will
752 // not be updated. We use this to only process internal edges once.
753 root_blossom_node_[subnode] = lca_index;
754
755 for (const EdgeIndex e : graph_[subnode]) {
756 Edge& edge = edges_[e];
757 const NodeIndex other_end = OtherEnd(edge, subnode);
758
759 // Skip edge that are already internal.
760 if (other_end == n) continue;
761
762 // This internal edge was already processed from its other end, so we
763 // can just skip it.
764 if (other_end == lca_index) {
765#ifndef NDEBUG
766 DCHECK_EQ(edge.slack, Slack(edge));
767#endif
768 continue;
769 }
770
771 // This is a new-internal edge that we didn't proccess yet.
772 //
773 // TODO(user): It would be nicer to not to have to read the memory of
774 // the other node at all. It might be possible once we store the
775 // parent edge instead of the parent node since then we will only need
776 // to know if this edges point to a new-internal node or not.
777 Node& mutable_other_node = nodes_[other_end];
778 if (mutable_other_node.is_internal) {
779 DCHECK(!plus_free_pq_.Contains(&edge));
780 if (plus_plus_pq_.Contains(&edge)) plus_plus_pq_.Remove(&edge);
781 edge.pseudo_slack += slack_adjust;
782 edge.pseudo_slack +=
783 mutable_other_node.IsMinus() ? tree_dual : -tree_dual;
784 continue;
785 }
786
787 // Replace the parent of any child of n by lca_index.
788 if (mutable_other_node.parent == n) {
789 mutable_other_node.parent = lca_index;
790 }
791
792 // Adjust when the edge used to be connected to a [-] node now that we
793 // attach it to a [+] node. Note that if the node was [+] then the
794 // non-internal incident edges slack and type do not change.
795 if (was_minus) {
796 edge.pseudo_slack += 2 * tree_dual;
797
798 // Add it to the correct PQ.
799 DCHECK(!plus_plus_pq_.Contains(&edge));
800 DCHECK(!plus_free_pq_.Contains(&edge));
801 if (mutable_other_node.IsPlus()) {
802 plus_plus_pq_.Add(&edge);
803 if (edge.pseudo_slack == 2 * tree_dual) {
804 primal_update_edge_queue_.push_back(e);
805 }
806 } else if (mutable_other_node.IsFree()) {
807 plus_free_pq_.Add(&edge);
808 if (edge.pseudo_slack == tree_dual) {
809 primal_update_edge_queue_.push_back(e);
810 }
811 }
812 }
813
814#ifndef NDEBUG
815 DCHECK_EQ(edge.slack, Slack(edge));
816#endif
817 }
818 }
819 }
820
821 DCHECK(backup_node.saved_blossom.empty());
822 backup_node.saved_blossom = std::move(lca.blossom);
823 lca.blossom = std::move(blossom);
824
825 VLOG(2) << "S result " << NodeDebugString(lca_index);
826}
827
828BlossomGraph::EdgeIndex BlossomGraph::FindTightExternalEdgeBetweenNodes(
831 DCHECK_EQ(tail, root_blossom_node_[tail]);
832 DCHECK_EQ(head, root_blossom_node_[head]);
833 for (const NodeIndex subnode : SubNodes(tail)) {
834 for (const EdgeIndex e : graph_[subnode]) {
835 const Edge& edge = edges_[e];
836 const NodeIndex other_end = OtherEnd(edge, subnode);
837 if (other_end == head && Slack(edge) == 0) {
838 return e;
839 }
840 }
841 }
842 return kNoEdgeIndex;
843}
844
846 ++num_expands_;
847 VLOG(2) << "Expand " << to_expand;
848
849 Node& node_to_expand = nodes_[to_expand];
850 DCHECK(node_to_expand.IsBlossom());
851 DCHECK(node_to_expand.IsMinus());
852 DCHECK_EQ(Dual(node_to_expand), 0);
853
854 const EdgeIndex match_edge_index =
855 FindTightExternalEdgeBetweenNodes(to_expand, node_to_expand.match);
856 const EdgeIndex parent_edge_index =
857 FindTightExternalEdgeBetweenNodes(to_expand, node_to_expand.parent);
858
859 // First, restore the saved fields.
860 Node& backup_node = nodes_[node_to_expand.blossom[1]];
861#ifndef NDEBUG
862 node_to_expand.dual = backup_node.saved_dual;
863#endif
864 node_to_expand.pseudo_dual = backup_node.saved_pseudo_dual;
865 std::vector<NodeIndex> blossom = std::move(node_to_expand.blossom);
866 node_to_expand.blossom = std::move(backup_node.saved_blossom);
867 backup_node.saved_blossom.clear();
868
869 // Restore the edges Head()/Tail().
870 for (const NodeIndex n : blossom) {
871 for (const NodeIndex subnode : SubNodes(n)) {
872 root_blossom_node_[subnode] = n;
873 }
874 }
875
876 // Now we try to find a 'blossom path' that will replace the blossom node in
877 // the alternating tree: the blossom's parent [+] node in the tree will be
878 // attached to a blossom subnode (the "path start"), the blossom's child in
879 // the tree will be attached to a blossom subnode (the "path end", which could
880 // be the same subnode or a different one), and, using the blossom cycle,
881 // we'll get a path with an odd number of blossom subnodes to connect the two
882 // (since the cycle is odd, one of the two paths will be odd too). The other
883 // subnodes of the blossom will then be made free nodes matched pairwise.
884 int blossom_path_start = -1;
885 int blossom_path_end = -1;
886 const NodeIndex start_node = OtherEndFromExternalNode(
887 edges_[parent_edge_index], node_to_expand.parent);
888 const NodeIndex end_node =
889 OtherEndFromExternalNode(edges_[match_edge_index], node_to_expand.match);
890 for (int i = 0; i < blossom.size(); ++i) {
891 if (blossom[i] == start_node) blossom_path_start = i;
892 if (blossom[i] == end_node) blossom_path_end = i;
893 }
894
895 // Split the cycle in two halves: nodes in [start..end] in path1, and
896 // nodes in [end..start] in path2. Note the inclusive intervals.
897 const std::vector<NodeIndex>& cycle = blossom;
898 std::vector<NodeIndex> path1;
899 std::vector<NodeIndex> path2;
900 {
901 const int end_offset =
902 (blossom_path_end + cycle.size() - blossom_path_start) % cycle.size();
903 for (int offset = 0; offset <= /*or equal*/ cycle.size(); ++offset) {
904 const NodeIndex node =
905 cycle[(blossom_path_start + offset) % cycle.size()];
906 if (offset <= end_offset) path1.push_back(node);
907 if (offset >= end_offset) path2.push_back(node);
908 }
909 }
910
911 // Reverse path2 to also make it go from start to end.
912 std::reverse(path2.begin(), path2.end());
913
914 // Swap if necessary so that path1 is the odd-length one.
915 if (path1.size() % 2 == 0) path1.swap(path2);
916
917 // Use better aliases than 'path1' and 'path2' in the code below.
918 std::vector<NodeIndex>& path_in_tree = path1;
919 const std::vector<NodeIndex>& free_pairs = path2;
920
921 // Strip path2 from the start and end, which aren't needed.
922 path2.erase(path2.begin());
923 path2.pop_back();
924
925 const NodeIndex blossom_matched_node = node_to_expand.match;
926 VLOG(2) << "Path ["
927 << absl::StrJoin(path_in_tree, ", ", absl::StreamFormatter())
928 << "] === " << blossom_matched_node;
929 VLOG(2) << "Pairs ["
930 << absl::StrJoin(free_pairs, ", ", absl::StreamFormatter()) << "]";
931
932 // Restore the path in the tree, note that we append the blossom_matched_node
933 // to simplify the code:
934 // <---- Blossom ---->
935 // [-] === [+] --- [-] === [+]
936 path_in_tree.push_back(blossom_matched_node);
937 CHECK_EQ(path_in_tree.size() % 2, 0);
938 const CostValue tree_dual = nodes_[node_to_expand.root].tree_dual_delta;
939 for (int i = 0; i < path_in_tree.size(); ++i) {
940 const NodeIndex n = path_in_tree[i];
941 const bool node_is_plus = i % 2;
942
943 // Update the parent.
944 if (i == 0) {
945 // This is the path start and its parent is either itself or the parent of
946 // to_expand if there was one.
947 DCHECK(node_to_expand.parent != to_expand || n == to_expand);
948 nodes_[n].parent = node_to_expand.parent;
949 } else {
950 nodes_[n].parent = path_in_tree[i - 1];
951 }
952
953 // Update the types and matches.
954 nodes_[n].root = node_to_expand.root;
955 nodes_[n].type = node_is_plus ? 1 : -1;
956 nodes_[n].match = path_in_tree[node_is_plus ? i - 1 : i + 1];
957
958 // Ignore the blossom_matched_node for the code below.
959 if (i + 1 == path_in_tree.size()) continue;
960
961 // Update the duals, depending on whether we have a new [+] or [-] node.
962 // Note that this is also needed for the 'root' blossom node (i=0), because
963 // we've restored its pseudo-dual from its old saved value above.
964 const CostValue adjust = node_is_plus ? -tree_dual : tree_dual;
965 nodes_[n].pseudo_dual += adjust;
966 for (const NodeIndex subnode : SubNodes(n)) {
967 for (const EdgeIndex e : graph_[subnode]) {
968 Edge& edge = edges_[e];
969 const NodeIndex other_end = OtherEnd(edge, subnode);
970 if (other_end == n) continue;
971
972 edge.pseudo_slack -= adjust;
973
974 // non-internal edges used to be attached to the [-] node_to_expand,
975 // so we adjust their dual.
976 if (other_end != to_expand && !nodes_[other_end].is_internal) {
977 edge.pseudo_slack += tree_dual;
978 } else {
979 // This was an internal edges. For the PQ code below to be correct, we
980 // wait for its other end to have been processed by this loop already.
981 // We detect that using the fact that the type of unprocessed internal
982 // node is still zero.
983 if (nodes_[other_end].type == 0) continue;
984 }
985
986 // Update edge queues.
987 if (node_is_plus) {
988 const Node& other_node = nodes_[other_end];
989 DCHECK(!plus_plus_pq_.Contains(&edge));
990 DCHECK(!plus_free_pq_.Contains(&edge));
991 if (other_node.IsPlus()) {
992 plus_plus_pq_.Add(&edge);
993 if (edge.pseudo_slack == 2 * tree_dual) {
994 primal_update_edge_queue_.push_back(e);
995 }
996 } else if (other_node.IsFree()) {
997 plus_free_pq_.Add(&edge);
998 if (edge.pseudo_slack == tree_dual) {
999 primal_update_edge_queue_.push_back(e);
1000 }
1001 }
1002 }
1003 }
1004 }
1005 }
1006
1007 // Update free nodes.
1008 for (const NodeIndex n : free_pairs) {
1009 nodes_[n].type = 0;
1010 nodes_[n].parent = n;
1011 nodes_[n].root = n;
1012
1013 // Update edges slack and priority queue for the adjacent edges.
1014 for (const NodeIndex subnode : SubNodes(n)) {
1015 for (const EdgeIndex e : graph_[subnode]) {
1016 Edge& edge = edges_[e];
1017 const NodeIndex other_end = OtherEnd(edge, subnode);
1018 if (other_end == n) continue;
1019
1020 // non-internal edges used to be attached to the [-] node_to_expand,
1021 // so we adjust their dual.
1022 if (other_end != to_expand && !nodes_[other_end].is_internal) {
1023 edge.pseudo_slack += tree_dual;
1024 }
1025
1026 // Update PQ. Note that since this was attached to a [-] node it cannot
1027 // be in any queue. We will also never process twice the same edge here.
1028 DCHECK(!plus_plus_pq_.Contains(&edge));
1029 DCHECK(!plus_free_pq_.Contains(&edge));
1030 if (nodes_[other_end].IsPlus()) {
1031 plus_free_pq_.Add(&edge);
1032 if (edge.pseudo_slack == tree_dual) {
1033 primal_update_edge_queue_.push_back(e);
1034 }
1035 }
1036 }
1037 }
1038 }
1039
1040 // Matches the free pair together.
1041 CHECK_EQ(free_pairs.size() % 2, 0);
1042 for (int i = 0; i < free_pairs.size(); i += 2) {
1043 nodes_[free_pairs[i]].match = free_pairs[i + 1];
1044 nodes_[free_pairs[i + 1]].match = free_pairs[i];
1045 }
1046
1047 // Mark all node as external. We do that last so we could easily detect old
1048 // internal edges that are now external.
1049 for (const NodeIndex n : blossom) {
1050 nodes_[n].is_internal = false;
1051 }
1052}
1053
1055 // Queue of blossoms to expand.
1056 std::vector<NodeIndex> queue;
1057 for (NodeIndex n(0); n < nodes_.size(); ++n) {
1058 Node& node = nodes_[n];
1059 if (node.IsInternal()) continue;
1060
1061 // When this is called, there should be no more trees.
1062 CHECK(node.IsFree());
1063
1064 if (node.IsBlossom()) queue.push_back(n);
1065 }
1066
1067 // TODO(user): remove duplication with expand?
1068 while (!queue.empty()) {
1069 const NodeIndex to_expand = queue.back();
1070 queue.pop_back();
1071
1072 Node& node_to_expand = nodes_[to_expand];
1073 DCHECK(node_to_expand.IsBlossom());
1074
1075 // Find the edge used to match to_expand with Match(to_expand).
1076 const EdgeIndex match_edge_index =
1077 FindTightExternalEdgeBetweenNodes(to_expand, node_to_expand.match);
1078
1079 // Restore the saved data.
1080 Node& backup_node = nodes_[node_to_expand.blossom[1]];
1081#ifndef NDEBUG
1082 node_to_expand.dual = backup_node.saved_dual;
1083#endif
1084 node_to_expand.pseudo_dual = backup_node.saved_pseudo_dual;
1085
1086 std::vector<NodeIndex> blossom = std::move(node_to_expand.blossom);
1087 node_to_expand.blossom = std::move(backup_node.saved_blossom);
1088 backup_node.saved_blossom.clear();
1089
1090 // Restore the edges Head()/Tail().
1091 for (const NodeIndex n : blossom) {
1092 for (const NodeIndex subnode : SubNodes(n)) {
1093 root_blossom_node_[subnode] = n;
1094 }
1095 }
1096
1097 // Find the index of matched_node in the blossom list.
1098 int internal_matched_index = -1;
1099 const NodeIndex matched_node = OtherEndFromExternalNode(
1100 edges_[match_edge_index], node_to_expand.match);
1101 const int size = blossom.size();
1102 for (int i = 0; i < size; ++i) {
1103 if (blossom[i] == matched_node) {
1104 internal_matched_index = i;
1105 break;
1106 }
1107 }
1108 CHECK_NE(internal_matched_index, -1);
1109
1110 // Amongst the node_to_expand.blossom nodes, internal_matched_index is
1111 // matched with external_matched_node and the other are matched together.
1112 std::vector<NodeIndex> free_pairs;
1113 for (int i = (internal_matched_index + 1) % size;
1114 i != internal_matched_index; i = (i + 1) % size) {
1115 free_pairs.push_back(blossom[i]);
1116 }
1117
1118 // Clear root/parent/type of all internal nodes.
1119 for (const NodeIndex to_clear : blossom) {
1120 nodes_[to_clear].type = 0;
1121 nodes_[to_clear].is_internal = false;
1122 nodes_[to_clear].parent = to_clear;
1123 nodes_[to_clear].root = to_clear;
1124 }
1125
1126 // Matches the internal node with external one.
1127 const NodeIndex external_matched_node = node_to_expand.match;
1128 const NodeIndex internal_matched_node = blossom[internal_matched_index];
1129 nodes_[internal_matched_node].match = external_matched_node;
1130 nodes_[external_matched_node].match = internal_matched_node;
1131
1132 // Matches the free pair together.
1133 CHECK_EQ(free_pairs.size() % 2, 0);
1134 for (int i = 0; i < free_pairs.size(); i += 2) {
1135 nodes_[free_pairs[i]].match = free_pairs[i + 1];
1136 nodes_[free_pairs[i + 1]].match = free_pairs[i];
1137 }
1138
1139 // Now that the expansion is done, add to the queue any sub-blossoms.
1140 for (const NodeIndex n : blossom) {
1141 if (nodes_[n].IsBlossom()) queue.push_back(n);
1142 }
1143 }
1144}
1145
1146const std::vector<NodeIndex>& BlossomGraph::SubNodes(NodeIndex n) {
1147 // This should be only called on an external node. However, in Shrink() we
1148 // mark the node as internal early, so we just make sure the node as no saved
1149 // blossom field here.
1150 DCHECK(nodes_[n].saved_blossom.empty());
1151
1152 // Expand all the inner nodes under the node n. This will not be n iff node is
1153 // is in fact a blossom.
1154 subnodes_ = {n};
1155 for (int i = 0; i < subnodes_.size(); ++i) {
1156 const Node& node = nodes_[subnodes_[i]];
1157
1158 // Since the first node in each list is always the node above, we just
1159 // skip it to avoid listing twice the nodes.
1160 if (!node.blossom.empty()) {
1161 subnodes_.insert(subnodes_.end(), node.blossom.begin() + 1,
1162 node.blossom.end());
1163 }
1164
1165 // We also need to recursively expand the sub-blossom nodes, which are (if
1166 // any) in the "saved_blossom" field of the first internal node of each
1167 // blossom. Since we iterate on all internal nodes here, we simply consult
1168 // the "saved_blossom" field of all subnodes, and it works the same.
1169 if (!node.saved_blossom.empty()) {
1170 subnodes_.insert(subnodes_.end(), node.saved_blossom.begin() + 1,
1171 node.saved_blossom.end());
1172 }
1173 }
1174 return subnodes_;
1175}
1176
1178 const Node& node = nodes_[n];
1179 if (node.is_internal) {
1180 return absl::StrCat("[I] #", n.value());
1181 }
1182 const std::string type = !NodeIsMatched(n) ? "[*]"
1183 : node.type == 1 ? "[+]"
1184 : node.type == -1 ? "[-]"
1185 : "[0]";
1186 return absl::StrCat(
1187 type, " #", n.value(), " dual: ", Dual(node).value(),
1188 " parent: ", node.parent.value(), " match: ", node.match.value(),
1189 " blossom: [", absl::StrJoin(node.blossom, ", ", absl::StreamFormatter()),
1190 "]");
1191}
1192
1193std::string BlossomGraph::EdgeDebugString(EdgeIndex e) const {
1194 const Edge& edge = edges_[e];
1195 if (nodes_[Tail(edge)].is_internal || nodes_[Head(edge)].is_internal) {
1196 return absl::StrCat(Tail(edge).value(), "<->", Head(edge).value(),
1197 " internal ");
1198 }
1199 return absl::StrCat(Tail(edge).value(), "<->", Head(edge).value(),
1200 " slack: ", Slack(edge).value());
1201}
1202
1203std::string BlossomGraph::DebugString() const {
1204 std::string result = "Graph:\n";
1205 for (NodeIndex n(0); n < nodes_.size(); ++n) {
1206 absl::StrAppend(&result, NodeDebugString(n), "\n");
1207 }
1208 for (EdgeIndex e(0); e < edges_.size(); ++e) {
1209 absl::StrAppend(&result, EdgeDebugString(e), "\n");
1210 }
1211 return result;
1212}
1213
1215#ifndef NDEBUG
1216 nodes_[n].dual += delta;
1217 for (const NodeIndex subnode : SubNodes(n)) {
1218 for (const EdgeIndex e : graph_[subnode]) {
1219 Edge& edge = edges_[e];
1220 const NodeIndex other_end = OtherEnd(edge, subnode);
1221 if (other_end == n) continue;
1222 edges_[e].slack -= delta;
1223 }
1224 }
1225#endif
1226}
1227
1229 const Node& tail_node = nodes_[Tail(edge)];
1230 const Node& head_node = nodes_[Head(edge)];
1231 CostValue slack = edge.pseudo_slack;
1232 if (Tail(edge) == Head(edge)) return slack; // Internal...
1233
1234 if (!tail_node.is_internal && !head_node.is_internal) {
1235 slack -= tail_node.type * nodes_[tail_node.root].tree_dual_delta +
1236 head_node.type * nodes_[head_node.root].tree_dual_delta;
1237 }
1238#ifndef NDEBUG
1239 DCHECK_EQ(slack, edge.slack) << tail_node.type << " " << head_node.type
1240 << " " << Tail(edge) << "<->" << Head(edge);
1241#endif
1242 return slack;
1243}
1244
1245// Returns the dual value of the given node (which might be a pseudo-node).
1247 const CostValue dual =
1248 node.pseudo_dual + node.type * nodes_[node.root].tree_dual_delta;
1249#ifndef NDEBUG
1250 DCHECK_EQ(dual, node.dual);
1251#endif
1252 return dual;
1253}
1254
1256 if (dual_objective_ == kint64max) return CostValue(kint64max);
1257 CHECK_EQ(dual_objective_ % 2, 0);
1258 return dual_objective_ / 2;
1259}
1260
1261void BlossomGraph::AddToDualObjective(CostValue delta) {
1262 CHECK_GE(delta, 0);
1263 dual_objective_ = CostValue(CapAdd(dual_objective_.value(), delta.value()));
1264}
1265
1267 VLOG(1) << "num_grows: " << num_grows_;
1268 VLOG(1) << "num_augments: " << num_augments_;
1269 VLOG(1) << "num_shrinks: " << num_shrinks_;
1270 VLOG(1) << "num_expands: " << num_expands_;
1271 VLOG(1) << "num_dual_updates: " << num_dual_updates_;
1272}
1273
1274} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define CHECK_GE(val1, val2)
Definition: base/logging.h:701
#define CHECK_GT(val1, val2)
Definition: base/logging.h:702
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define CHECK_NE(val1, val2)
Definition: base/logging.h:698
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
#define VLOG(verboselevel)
Definition: base/logging.h:978
void resize(size_type new_size)
void push_back(const value_type &x)
ABSL_MUST_USE_RESULT bool Initialize()
void Grow(EdgeIndex e, NodeIndex tail, NodeIndex head)
CostValue ComputeMaxCommonTreeDualDeltaAndResetPrimalEdgeQueue()
bool DebugEdgeIsTightAndExternal(const Edge &edge) const
static const CostValue kMaxCostValue
NodeIndex Match(NodeIndex n) const
void AddEdge(NodeIndex tail, NodeIndex head, CostValue cost)
bool NodeIsMatched(NodeIndex n) const
CostValue Slack(const Edge &edge) const
std::string NodeDebugString(NodeIndex n) const
std::string EdgeDebugString(EdgeIndex e) const
CostValue Dual(const Node &node) const
static const EdgeIndex kNoEdgeIndex
static const NodeIndex kNoNodeIndex
void Expand(NodeIndex to_expand)
void UpdateAllTrees(CostValue delta)
void DebugUpdateNodeDual(NodeIndex n, CostValue delta)
void AddEdgeWithCost(int tail, int head, int64 cost)
int64 value
static const int64 kint64max
int64_t int64
const int FATAL
Definition: log_severity.h:32
const bool DEBUG_MODE
Definition: macros.h:24
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapAdd(int64 x, int64 y)
int index
Definition: pack.cc:508
int64 delta
Definition: resource.cc:1684
int64 tail
int64 cost
int64 head
NodeIndex OtherEnd(NodeIndex n) const