OR-Tools  8.1
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 
19 namespace operations_research {
20 
21 void 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 
111 const 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 
538 void 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 
547 void 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 
639 int 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 
650 void 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 
828 BlossomGraph::EdgeIndex BlossomGraph::FindTightExternalEdgeBetweenNodes(
830  DCHECK_NE(tail, head);
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 
1146 const 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 
1193 std::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 
1203 std::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 
1228 CostValue BlossomGraph::Slack(const Edge& edge) const {
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).
1246 CostValue BlossomGraph::Dual(const Node& node) const {
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 
1261 void 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
tail
int64 tail
Definition: routing_flow.cc:127
operations_research::BlossomGraph::Node::saved_pseudo_dual
CostValue saved_pseudo_dual
Definition: perfect_matching.h:252
operations_research::MinCostPerfectMatching::Solve
ABSL_MUST_USE_RESULT Status Solve()
Definition: perfect_matching.cc:39
min
int64 min
Definition: alldiff_cst.cc:138
DCHECK_LT
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
VLOG
#define VLOG(verboselevel)
Definition: base/logging.h:978
absl::StrongVector::push_back
void push_back(const value_type &x)
Definition: strong_vector.h:158
operations_research::BlossomGraph::Slack
CostValue Slack(const Edge &edge) const
Definition: perfect_matching.cc:1228
operations_research::BlossomGraph::NodeDebugString
std::string NodeDebugString(NodeIndex n) const
Definition: perfect_matching.cc:1177
operations_research::BlossomGraph::Node::IsMinus
bool IsMinus() const
Definition: perfect_matching.h:195
max
int64 max
Definition: alldiff_cst.cc:139
LOG
#define LOG(severity)
Definition: base/logging.h:420
operations_research::BlossomGraph::DebugCheckNoPossiblePrimalUpdates
void DebugCheckNoPossiblePrimalUpdates()
Definition: perfect_matching.cc:347
operations_research::NodeIndex
int32 NodeIndex
Definition: ebert_graph.h:192
FATAL
const int FATAL
Definition: log_severity.h:32
operations_research::BlossomGraph::Edge::OtherEnd
NodeIndex OtherEnd(NodeIndex n) const
Definition: perfect_matching.h:268
operations_research::BlossomGraph::Node::IsBlossom
bool IsBlossom() const
Definition: perfect_matching.h:200
CHECK_GE
#define CHECK_GE(val1, val2)
Definition: base/logging.h:701
operations_research::BlossomGraph::EdgeDebugString
std::string EdgeDebugString(EdgeIndex e) const
Definition: perfect_matching.cc:1193
operations_research::BlossomGraph::DebugUpdateNodeDual
void DebugUpdateNodeDual(NodeIndex n, CostValue delta)
Definition: perfect_matching.cc:1214
value
int64 value
Definition: demon_profiler.cc:43
CHECK_GT
#define CHECK_GT(val1, val2)
Definition: base/logging.h:702
operations_research::BlossomGraph::Node::saved_blossom
std::vector< NodeIndex > saved_blossom
Definition: perfect_matching.h:253
operations_research::BlossomGraph::BlossomGraph
BlossomGraph(int num_nodes)
Definition: perfect_matching.cc:116
operations_research::BlossomGraph::Node::parent
NodeIndex parent
Definition: perfect_matching.h:214
saturated_arithmetic.h
operations_research
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
Definition: dense_doubly_linked_list.h:21
operations_research::BlossomGraph::Expand
void Expand(NodeIndex to_expand)
Definition: perfect_matching.cc:845
operations_research::BlossomGraph::UpdateAllTrees
void UpdateAllTrees(CostValue delta)
Definition: perfect_matching.cc:309
operations_research::BlossomGraph::Node::type
int type
Definition: perfect_matching.h:207
operations_research::BlossomGraph::Edge
Definition: perfect_matching.h:257
int64
int64_t int64
Definition: integral_types.h:34
operations_research::BlossomGraph::ComputeMaxCommonTreeDualDeltaAndResetPrimalEdgeQueue
CostValue ComputeMaxCommonTreeDualDeltaAndResetPrimalEdgeQueue()
Definition: perfect_matching.cc:259
operations_research::MinCostPerfectMatching::Reset
void Reset(int num_nodes)
Definition: perfect_matching.cc:21
operations_research::BlossomGraph::PrimalUpdates
void PrimalUpdates()
Definition: perfect_matching.cc:381
index
int index
Definition: pack.cc:508
operations_research::BlossomGraph::DebugEdgeIsTightAndExternal
bool DebugEdgeIsTightAndExternal(const Edge &edge) const
Definition: perfect_matching.cc:463
operations_research::BlossomGraph::Node::match
NodeIndex match
Definition: perfect_matching.h:218
operations_research::BlossomGraph::Node::dual
CostValue dual
Definition: perfect_matching.h:235
DCHECK_NE
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
cost
int64 cost
Definition: routing_flow.cc:130
operations_research::BlossomGraph::Node::pseudo_dual
CostValue pseudo_dual
Definition: perfect_matching.h:231
DEBUG_MODE
const bool DEBUG_MODE
Definition: macros.h:24
operations_research::MinCostPerfectMatching::AddEdgeWithCost
void AddEdgeWithCost(int tail, int head, int64 cost)
Definition: perfect_matching.cc:27
operations_research::sat::INFEASIBLE
@ INFEASIBLE
Definition: cp_model.pb.h:226
operations_research::BlossomGraph::Node
Definition: perfect_matching.h:181
operations_research::BlossomGraph::Node::IsFree
bool IsFree() const
Definition: perfect_matching.h:193
operations_research::BlossomGraph::Shrink
void Shrink(EdgeIndex e)
Definition: perfect_matching.cc:650
operations_research::MinCostPerfectMatching::Status
Status
Definition: perfect_matching.h:81
operations_research::CapAdd
int64 CapAdd(int64 x, int64 y)
Definition: saturated_arithmetic.h:124
operations_research::CostValue
int64 CostValue
Definition: ebert_graph.h:203
operations_research::BlossomGraph::ExpandAllBlossoms
void ExpandAllBlossoms()
Definition: perfect_matching.cc:1054
operations_research::BlossomGraph::DualObjective
CostValue DualObjective() const
Definition: perfect_matching.cc:1255
operations_research::BlossomGraph::kNoEdgeIndex
static const EdgeIndex kNoEdgeIndex
Definition: perfect_matching.h:176
CHECK_EQ
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
operations_research::BlossomGraph::Initialize
ABSL_MUST_USE_RESULT bool Initialize()
Definition: perfect_matching.cc:143
operations_research::BlossomGraph::Edge::slack
CostValue slack
Definition: perfect_matching.h:286
operations_research::sat::OPTIMAL
@ OPTIMAL
Definition: cp_model.pb.h:227
operations_research::BlossomGraph::DebugDualsAreFeasible
bool DebugDualsAreFeasible() const
Definition: perfect_matching.cc:450
DCHECK
#define DCHECK(condition)
Definition: base/logging.h:884
operations_research::BlossomGraph::Node::root
NodeIndex root
Definition: perfect_matching.h:222
operations_research::BlossomGraph::Edge::head
NodeIndex head
Definition: perfect_matching.h:295
operations_research::BlossomGraph::Edge::pseudo_slack
CostValue pseudo_slack
Definition: perfect_matching.h:282
DCHECK_GE
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
operations_research::BlossomGraph::Node::IsPlus
bool IsPlus() const
Definition: perfect_matching.h:194
operations_research::BlossomGraph::Match
NodeIndex Match(NodeIndex n) const
Definition: perfect_matching.cc:336
operations_research::BlossomGraph::DebugString
std::string DebugString() const
Definition: perfect_matching.cc:1203
operations_research::BlossomGraph::AddEdge
void AddEdge(NodeIndex tail, NodeIndex head, CostValue cost)
Definition: perfect_matching.cc:126
operations_research::BlossomGraph::kMaxCostValue
static const CostValue kMaxCostValue
Definition: perfect_matching.h:177
DCHECK_EQ
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
operations_research::BlossomGraph::Augment
void Augment(EdgeIndex e)
Definition: perfect_matching.cc:547
delta
int64 delta
Definition: resource.cc:1684
absl::StrongVector::resize
void resize(size_type new_size)
Definition: strong_vector.h:150
operations_research::BlossomGraph::Node::IsInternal
bool IsInternal() const
Definition: perfect_matching.h:189
CHECK_NE
#define CHECK_NE(val1, val2)
Definition: base/logging.h:698
operations_research::BlossomGraph::NodeIsMatched
bool NodeIsMatched(NodeIndex n) const
Definition: perfect_matching.cc:329
operations_research::BlossomGraph::Dual
CostValue Dual(const Node &node) const
Definition: perfect_matching.cc:1246
head
int64 head
Definition: routing_flow.cc:128
operations_research::BlossomGraph::Grow
void Grow(EdgeIndex e, NodeIndex tail, NodeIndex head)
Definition: perfect_matching.cc:470
operations_research::BlossomGraph::kNoNodeIndex
static const NodeIndex kNoNodeIndex
Definition: perfect_matching.h:175
DCHECK_LE
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
operations_research::BlossomGraph::Edge::tail
NodeIndex tail
Definition: perfect_matching.h:294
CHECK
#define CHECK(condition)
Definition: base/logging.h:495
operations_research::BlossomGraph::Node::saved_dual
CostValue saved_dual
Definition: perfect_matching.h:250
operations_research::BlossomGraph::DisplayStats
void DisplayStats() const
Definition: perfect_matching.cc:1266
kint64max
static const int64 kint64max
Definition: integral_types.h:62
perfect_matching.h
operations_research::BlossomGraph::Node::blossom
std::vector< NodeIndex > blossom
Definition: perfect_matching.h:241
operations_research::BlossomGraph::Node::is_internal
bool is_internal
Definition: perfect_matching.h:210