OR-Tools  8.2
linear_assignment.h
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
14//
15// An implementation of a cost-scaling push-relabel algorithm for the
16// assignment problem (minimum-cost perfect bipartite matching), from
17// the paper of Goldberg and Kennedy (1995).
18//
19//
20// This implementation finds the minimum-cost perfect assignment in
21// the given graph with integral edge weights set through the
22// SetArcCost method.
23//
24// The running time is O(n*m*log(nC)) where n is the number of nodes,
25// m is the number of edges, and C is the largest magnitude of an edge cost.
26// In principle it can be worse than the Hungarian algorithm but we don't know
27// of any class of problems where that actually happens. An additional sqrt(n)
28// factor could be shaved off the running time bound using the technique
29// described in http://dx.doi.org/10.1137/S0895480194281185
30// (see also http://theory.stanford.edu/~robert/papers/glob_upd.ps).
31//
32//
33// Example usage:
34//
35// #include "ortools/graph/graph.h"
36// #include "ortools/graph/linear_assignment.h"
37//
38// // Choose a graph implementation (we recommend StaticGraph<>).
39// typedef util::StaticGraph<> Graph;
40//
41// // Define a num_nodes / 2 by num_nodes / 2 assignment problem:
42// const int num_nodes = ...
43// const int num_arcs = ...
44// const int num_left_nodes = num_nodes / 2;
45// Graph graph(num_nodes, num_arcs);
46// std::vector<operations_research::CostValue> arc_costs(num_arcs);
47// for (int arc = 0; arc < num_arcs; ++arc) {
48// const int arc_tail = ... // must be in [0, num_left_nodes)
49// const int arc_head = ... // must be in [num_left_nodes, num_nodes)
50// graph.AddArc(arc_tail, arc_head);
51// arc_costs[arc] = ...
52// }
53//
54// // Build the StaticGraph. You can skip this step by using a ListGraph<>
55// // instead, but then the ComputeAssignment() below will be slower. It is
56// // okay if your graph is small and performance is not critical though.
57// {
58// std::vector<Graph::ArcIndex> arc_permutation;
59// graph.Build(&arc_permutation);
60// util::Permute(arc_permutation, &arc_costs);
61// }
62//
63// // Construct the LinearSumAssignment.
64// ::operations_research::LinearSumAssignment<Graph> a(graph, num_left_nodes);
65// for (int arc = 0; arc < num_arcs; ++arc) {
66// // You can also replace 'arc_costs[arc]' by something like
67// // ComputeArcCost(permutation.empty() ? arc : permutation[arc])
68// // if you don't want to store the costs in arc_costs to save memory.
69// a.SetArcCost(arc, arc_costs[arc]);
70// }
71//
72// // Compute the optimum assignment.
73// bool success = a.ComputeAssignment();
74// // Retrieve the cost of the optimum assignment.
75// operations_research::CostValue optimum_cost = a.GetCost();
76// // Retrieve the node-node correspondence of the optimum assignment and the
77// // cost of each node pairing.
78// for (int left_node = 0; left_node < num_left_nodes; ++left_node) {
79// const int right_node = a.GetMate(left_node);
80// operations_research::CostValue node_pair_cost =
81// a.GetAssignmentCost(left_node);
82// ...
83// }
84//
85// In the following, we consider a bipartite graph
86// G = (V = X union Y, E subset XxY),
87// where V denotes the set of nodes (vertices) in the graph, E denotes
88// the set of arcs (edges), n = |V| denotes the number of nodes in the
89// graph, and m = |E| denotes the number of arcs in the graph.
90//
91// The set of nodes is divided into two parts, X and Y, and every arc
92// must go between a node of X and a node of Y. With each arc is
93// associated a cost c(v, w). A matching M is a subset of E with the
94// property that no two arcs in M have a head or tail node in common,
95// and a perfect matching is a matching that touches every node in the
96// graph. The cost of a matching M is the sum of the costs of all the
97// arcs in M.
98//
99// The assignment problem is to find a perfect matching of minimum
100// cost in the given bipartite graph. The present algorithm reduces
101// the assignment problem to an instance of the minimum-cost flow
102// problem and takes advantage of special properties of the resulting
103// minimum-cost flow problem to solve it efficiently using a
104// push-relabel method. For more information about minimum-cost flow
105// see google3/ortools/graph/min_cost_flow.h
106//
107// The method used here is the cost-scaling approach for the
108// minimum-cost circulation problem as described in [Goldberg and
109// Tarjan] with some technical modifications:
110// 1. For efficiency, we solve a transportation problem instead of
111// minimum-cost circulation. We might revisit this decision if it
112// is important to handle problems in which no perfect matching
113// exists.
114// 2. We use a modified "asymmetric" notion of epsilon-optimality in
115// which left-to-right residual arcs are required to have reduced
116// cost bounded below by zero and right-to-left residual arcs are
117// required to have reduced cost bounded below by -epsilon. For
118// each residual arc direction, the reduced-cost threshold for
119// admissibility is epsilon/2 above the threshold for epsilon
120// optimality.
121// 3. We do not limit the applicability of the relabeling operation to
122// nodes with excess. Instead we use the double-push operation
123// (discussed in the Goldberg and Kennedy CSA paper and Kennedy's
124// thesis) which relabels right-side nodes just *after* they have
125// been discharged.
126// The above differences are explained in detail in [Kennedy's thesis]
127// and explained not quite as cleanly in [Goldberg and Kennedy's CSA
128// paper]. But note that the thesis explanation uses a value of
129// epsilon that's double what we use here.
130//
131// Some definitions:
132// Active: A node is called active when it has excess. It is
133// eligible to be pushed from. In this implementation, every active
134// node is on the left side of the graph where prices are determined
135// implicitly, so no left-side relabeling is necessary before
136// pushing from an active node. We do, however, need to compute
137// the implications for price changes on the affected right-side
138// nodes.
139// Admissible: A residual arc (one that can carry more flow) is
140// called admissible when its reduced cost is small enough. We can
141// push additional flow along such an arc without violating
142// epsilon-optimality. In the case of a left-to-right residual
143// arc, the reduced cost must be at most epsilon/2. In the case of
144// a right-to-left residual arc, the reduced cost must be at
145// most -epsilon/2. The careful reader will note that these thresholds
146// are not used explicitly anywhere in this implementation, and
147// the reason is the implicit pricing of left-side nodes.
148// Reduced cost: Essentially an arc's reduced cost is its
149// complementary slackness. In push-relabel algorithms this is
150// c_p(v, w) = p(v) + c(v, w) - p(w),
151// where p() is the node price function and c(v, w) is the cost of
152// the arc from v to w. See min_cost_flow.h for more details.
153// Partial reduced cost: We maintain prices implicitly for left-side
154// nodes in this implementation, so instead of reduced costs we
155// work with partial reduced costs, defined as
156// c'_p(v, w) = c(v, w) - p(w).
157//
158// We check at initialization time for the possibility of arithmetic
159// overflow and warn if the given costs are too large. In many cases
160// the bound we use to trigger the warning is pessimistic so the given
161// problem can often be solved even if we warn that overflow is
162// possible.
163//
164// We don't use the interface from
165// operations_research/algorithms/hungarian.h because we want to be
166// able to express sparse problems efficiently.
167//
168// When asked to solve the given assignment problem we return a
169// boolean to indicate whether the given problem was feasible.
170//
171// References:
172// [ Goldberg and Kennedy's CSA paper ] A. V. Goldberg and R. Kennedy,
173// "An Efficient Cost Scaling Algorithm for the Assignment Problem."
174// Mathematical Programming, Vol. 71, pages 153-178, December 1995.
175//
176// [ Goldberg and Tarjan ] A. V. Goldberg and R. E. Tarjan, "Finding
177// Minimum-Cost Circulations by Successive Approximation." Mathematics
178// of Operations Research, Vol. 15, No. 3, pages 430-466, August 1990.
179//
180// [ Kennedy's thesis ] J. R. Kennedy, Jr., "Solving Unweighted and
181// Weighted Bipartite Matching Problems in Theory and Practice."
182// Stanford University Doctoral Dissertation, Department of Computer
183// Science, 1995.
184//
185// [ Burkard et al. ] R. Burkard, M. Dell'Amico, S. Martello, "Assignment
186// Problems", SIAM, 2009, ISBN: 978-0898716634,
187// http://www.amazon.com/dp/0898716632/
188//
189// [ Ahuja et al. ] R. K. Ahuja, T. L. Magnanti, J. B. Orlin, "Network Flows:
190// Theory, Algorithms, and Applications," Prentice Hall, 1993,
191// ISBN: 978-0136175490, http://www.amazon.com/dp/013617549X.
192//
193// Keywords: linear sum assignment problem, Hungarian method, Goldberg, Kennedy.
194
195#ifndef OR_TOOLS_GRAPH_LINEAR_ASSIGNMENT_H_
196#define OR_TOOLS_GRAPH_LINEAR_ASSIGNMENT_H_
197
198#include <algorithm>
199#include <cstdlib>
200#include <deque>
201#include <limits>
202#include <memory>
203#include <string>
204#include <utility>
205#include <vector>
206
207#include "absl/strings/str_format.h"
210#include "ortools/base/logging.h"
211#include "ortools/base/macros.h"
214#include "ortools/util/zvector.h"
215
216#ifndef SWIG
217ABSL_DECLARE_FLAG(int64, assignment_alpha);
218ABSL_DECLARE_FLAG(int, assignment_progress_logging_period);
219ABSL_DECLARE_FLAG(bool, assignment_stack_order);
220#endif
221
222namespace operations_research {
223
224// This class does not take ownership of its underlying graph.
225template <typename GraphType>
227 public:
230
231 // Constructor for the case in which we will build the graph
232 // incrementally as we discover arc costs, as might be done with any
233 // of the dynamic graph representations such as StarGraph or ForwardStarGraph.
234 LinearSumAssignment(const GraphType& graph, NodeIndex num_left_nodes);
235
236 // Constructor for the case in which the underlying graph cannot be
237 // built until after all the arc costs are known, as is the case
238 // with ForwardStarStaticGraph. In this case, the graph is passed to
239 // us later via the SetGraph() method, below.
240 LinearSumAssignment(NodeIndex num_left_nodes, ArcIndex num_arcs);
241
243
244 // Sets the graph used by the LinearSumAssignment instance, for use
245 // when the graph layout can be determined only after arc costs are
246 // set. This happens, for example, when we use a ForwardStarStaticGraph.
247 void SetGraph(const GraphType* graph) {
248 DCHECK(graph_ == nullptr);
249 graph_ = graph;
250 }
251
252 // Sets the cost-scaling divisor, i.e., the amount by which we
253 // divide the scaling parameter on each iteration.
254 void SetCostScalingDivisor(CostValue factor) { alpha_ = factor; }
255
256 // Returns a permutation cycle handler that can be passed to the
257 // TransformToForwardStaticGraph method so that arc costs get
258 // permuted along with arcs themselves.
259 //
260 // Passes ownership of the cycle handler to the caller.
261 //
264
265 // Optimizes the layout of the graph for the access pattern our
266 // implementation will use.
267 //
268 // REQUIRES for LinearSumAssignment template instantiation if a call
269 // to the OptimizeGraphLayout() method is compiled: GraphType is a
270 // dynamic graph, i.e., one that implements the
271 // GroupForwardArcsByFunctor() member template method.
272 //
273 // If analogous optimization is needed for LinearSumAssignment
274 // instances based on static graphs, the graph layout should be
275 // constructed such that each node's outgoing arcs are sorted by
276 // head node index before the
277 // LinearSumAssignment<GraphType>::SetGraph() method is called.
278 void OptimizeGraphLayout(GraphType* graph);
279
280 // Allows tests, iterators, etc., to inspect our underlying graph.
281 inline const GraphType& Graph() const { return *graph_; }
282
283 // These handy member functions make the code more compact, and we
284 // expose them to clients so that client code that doesn't have
285 // direct access to the graph can learn about the optimum assignment
286 // once it is computed.
287 inline NodeIndex Head(ArcIndex arc) const { return graph_->Head(arc); }
288
289 // Returns the original arc cost for use by a client that's
290 // iterating over the optimum assignment.
292 DCHECK_EQ(0, scaled_arc_cost_[arc] % cost_scaling_factor_);
293 return scaled_arc_cost_[arc] / cost_scaling_factor_;
294 }
295
296 // Sets the cost of an arc already present in the given graph.
298
299 // Completes initialization after the problem is fully specified.
300 // Returns true if we successfully prove that arithmetic
301 // calculations are guaranteed not to overflow. ComputeAssignment()
302 // calls this method itself, so only clients that care about
303 // obtaining a warning about the possibility of arithmetic precision
304 // problems need to call this method explicitly.
305 //
306 // Separate from ComputeAssignment() for white-box testing and for
307 // clients that need to react to the possibility that arithmetic
308 // overflow is not ruled out.
309 //
310 // FinalizeSetup() is idempotent.
311 bool FinalizeSetup();
312
313 // Computes the optimum assignment. Returns true on success. Return
314 // value of false implies the given problem is infeasible.
315 bool ComputeAssignment();
316
317 // Returns the cost of the minimum-cost perfect matching.
318 // Precondition: success_ == true, signifying that we computed the
319 // optimum assignment for a feasible problem.
320 CostValue GetCost() const;
321
322 // Returns the total number of nodes in the given problem.
324 if (graph_ == nullptr) {
325 // Return a guess that must be true if ultimately we are given a
326 // feasible problem to solve.
327 return 2 * NumLeftNodes();
328 } else {
329 return graph_->num_nodes();
330 }
331 }
332
333 // Returns the number of nodes on the left side of the given
334 // problem.
335 NodeIndex NumLeftNodes() const { return num_left_nodes_; }
336
337 // Returns the arc through which the given node is matched.
338 inline ArcIndex GetAssignmentArc(NodeIndex left_node) const {
339 DCHECK_LT(left_node, num_left_nodes_);
340 return matched_arc_[left_node];
341 }
342
343 // Returns the cost of the assignment arc incident to the given
344 // node.
346 return ArcCost(GetAssignmentArc(node));
347 }
348
349 // Returns the node to which the given node is matched.
350 inline NodeIndex GetMate(NodeIndex left_node) const {
351 DCHECK_LT(left_node, num_left_nodes_);
352 ArcIndex matching_arc = GetAssignmentArc(left_node);
353 DCHECK_NE(GraphType::kNilArc, matching_arc);
354 return Head(matching_arc);
355 }
356
357 std::string StatsString() const { return total_stats_.StatsString(); }
358
360 public:
361 BipartiteLeftNodeIterator(const GraphType& graph, NodeIndex num_left_nodes)
362 : num_left_nodes_(num_left_nodes), node_iterator_(0) {}
363
365 : num_left_nodes_(assignment.NumLeftNodes()), node_iterator_(0) {}
366
367 NodeIndex Index() const { return node_iterator_; }
368
369 bool Ok() const { return node_iterator_ < num_left_nodes_; }
370
371 void Next() { ++node_iterator_; }
372
373 private:
374 const NodeIndex num_left_nodes_;
375 typename GraphType::NodeIndex node_iterator_;
376 };
377
378 private:
379 struct Stats {
380 Stats() : pushes_(0), double_pushes_(0), relabelings_(0), refinements_(0) {}
381 void Clear() {
382 pushes_ = 0;
383 double_pushes_ = 0;
384 relabelings_ = 0;
385 refinements_ = 0;
386 }
387 void Add(const Stats& that) {
388 pushes_ += that.pushes_;
389 double_pushes_ += that.double_pushes_;
390 relabelings_ += that.relabelings_;
391 refinements_ += that.refinements_;
392 }
393 std::string StatsString() const {
394 return absl::StrFormat(
395 "%d refinements; %d relabelings; "
396 "%d double pushes; %d pushes",
397 refinements_, relabelings_, double_pushes_, pushes_);
398 }
399 int64 pushes_;
400 int64 double_pushes_;
401 int64 relabelings_;
402 int64 refinements_;
403 };
404
405#ifndef SWIG
406 class ActiveNodeContainerInterface {
407 public:
408 virtual ~ActiveNodeContainerInterface() {}
409 virtual bool Empty() const = 0;
410 virtual void Add(NodeIndex node) = 0;
411 virtual NodeIndex Get() = 0;
412 };
413
414 class ActiveNodeStack : public ActiveNodeContainerInterface {
415 public:
416 ~ActiveNodeStack() override {}
417
418 bool Empty() const override { return v_.empty(); }
419
420 void Add(NodeIndex node) override { v_.push_back(node); }
421
422 NodeIndex Get() override {
423 DCHECK(!Empty());
424 NodeIndex result = v_.back();
425 v_.pop_back();
426 return result;
427 }
428
429 private:
430 std::vector<NodeIndex> v_;
431 };
432
433 class ActiveNodeQueue : public ActiveNodeContainerInterface {
434 public:
435 ~ActiveNodeQueue() override {}
436
437 bool Empty() const override { return q_.empty(); }
438
439 void Add(NodeIndex node) override { q_.push_front(node); }
440
441 NodeIndex Get() override {
442 DCHECK(!Empty());
443 NodeIndex result = q_.back();
444 q_.pop_back();
445 return result;
446 }
447
448 private:
449 std::deque<NodeIndex> q_;
450 };
451#endif
452
453 // Type definition for a pair
454 // (arc index, reduced cost gap)
455 // giving the arc along which we will push from a given left-side
456 // node and the gap between that arc's partial reduced cost and the
457 // reduced cost of the next-best (necessarily residual) arc out of
458 // the node. This information helps us efficiently relabel
459 // right-side nodes during DoublePush operations.
460 typedef std::pair<ArcIndex, CostValue> ImplicitPriceSummary;
461
462 // Returns true if and only if the current pseudoflow is
463 // epsilon-optimal. To be used in a DCHECK.
464 bool EpsilonOptimal() const;
465
466 // Checks that all nodes are matched.
467 // To be used in a DCHECK.
468 bool AllMatched() const;
469
470 // Calculates the implicit price of the given node.
471 // Only for debugging, for use in EpsilonOptimal().
472 inline CostValue ImplicitPrice(NodeIndex left_node) const;
473
474 // For use by DoublePush()
475 inline ImplicitPriceSummary BestArcAndGap(NodeIndex left_node) const;
476
477 // Accumulates stats between iterations and reports them if the
478 // verbosity level is high enough.
479 void ReportAndAccumulateStats();
480
481 // Utility function to compute the next error parameter value. This
482 // is used to ensure that the same sequence of error parameter
483 // values is used for computation of price bounds as is used for
484 // computing the optimum assignment.
485 CostValue NewEpsilon(CostValue current_epsilon) const;
486
487 // Advances internal state to prepare for the next scaling
488 // iteration. Returns false if infeasibility is detected, true
489 // otherwise.
490 bool UpdateEpsilon();
491
492 // Indicates whether the given left_node has positive excess. Called
493 // only for nodes on the left side.
494 inline bool IsActive(NodeIndex left_node) const;
495
496 // Indicates whether the given node has nonzero excess. The idea
497 // here is the same as the IsActive method above, but that method
498 // contains a safety DCHECK() that its argument is a left-side node,
499 // while this method is usable for any node.
500 // To be used in a DCHECK.
501 inline bool IsActiveForDebugging(NodeIndex node) const;
502
503 // Performs the push/relabel work for one scaling iteration.
504 bool Refine();
505
506 // Puts all left-side nodes in the active set in preparation for the
507 // first scaling iteration.
508 void InitializeActiveNodeContainer();
509
510 // Saturates all negative-reduced-cost arcs at the beginning of each
511 // scaling iteration. Note that according to the asymmetric
512 // definition of admissibility, this action is different from
513 // saturating all admissible arcs (which we never do). All negative
514 // arcs are admissible, but not all admissible arcs are negative. It
515 // is alwsys enough to saturate only the negative ones.
516 void SaturateNegativeArcs();
517
518 // Performs an optimized sequence of pushing a unit of excess out of
519 // the left-side node v and back to another left-side node if no
520 // deficit is cancelled with the first push.
521 bool DoublePush(NodeIndex source);
522
523 // Returns the partial reduced cost of the given arc.
524 inline CostValue PartialReducedCost(ArcIndex arc) const {
525 return scaled_arc_cost_[arc] - price_[Head(arc)];
526 }
527
528 // The graph underlying the problem definition we are given. Not
529 // owned by *this.
530 const GraphType* graph_;
531
532 // The number of nodes on the left side of the graph we are given.
533 NodeIndex num_left_nodes_;
534
535 // A flag indicating, after FinalizeSetup() has run, whether the
536 // arc-incidence precondition required by BestArcAndGap() is
537 // satisfied by every left-side node. If not, the problem is
538 // infeasible.
539 bool incidence_precondition_satisfied_;
540
541 // A flag indicating that an optimal perfect matching has been computed.
542 bool success_;
543
544 // The value by which we multiply all the arc costs we are given in
545 // order to be able to use integer arithmetic in all our
546 // computations. In order to establish optimality of the final
547 // matching we compute, we need that
548 // (cost_scaling_factor_ / kMinEpsilon) > graph_->num_nodes().
549 const CostValue cost_scaling_factor_;
550
551 // Scaling divisor.
552 CostValue alpha_;
553
554 // Minimum value of epsilon. When a flow is epsilon-optimal for
555 // epsilon == kMinEpsilon, the flow is optimal.
556 static const CostValue kMinEpsilon;
557
558 // Current value of epsilon, the cost scaling parameter.
559 CostValue epsilon_;
560
561 // The following two data members, price_lower_bound_ and
562 // slack_relabeling_price_, have to do with bounds on the amount by
563 // which node prices can change during execution of the algorithm.
564 // We need some detailed discussion of this topic because we violate
565 // several simplifying assumptions typically made in the theoretical
566 // literature. In particular, we use integer arithmetic, we use a
567 // reduction to the transportation problem rather than min-cost
568 // circulation, we provide detection of infeasible problems rather
569 // than assume feasibility, we detect when our computations might
570 // exceed the range of representable cost values, and we use the
571 // double-push heuristic which relabels nodes that do not have
572 // excess.
573 //
574 // In the following discussion, we prove the following propositions:
575 // Proposition 1. [Fidelity of arithmetic precision guarantee] If
576 // FinalizeSetup() returns true, no arithmetic
577 // overflow occurs during ComputeAssignment().
578 // Proposition 2. [Fidelity of feasibility detection] If no
579 // arithmetic overflow occurs during
580 // ComputeAssignment(), the return value of
581 // ComputeAssignment() faithfully indicates whether
582 // the given problem is feasible.
583 //
584 // We begin with some general discussion.
585 //
586 // The ideas used to prove our two propositions are essentially
587 // those that appear in [Goldberg and Tarjan], but several details
588 // are different: [Goldberg and Tarjan] assumes a feasible problem,
589 // uses a symmetric notion of epsilon-optimality, considers only
590 // nodes with excess eligible for relabeling, and does not treat the
591 // question of arithmetic overflow. This implementation, on the
592 // other hand, detects and reports infeasible problems, uses
593 // asymmetric epsilon-optimality, relabels nodes with no excess in
594 // the course of the double-push operation, and gives a reasonably
595 // tight guarantee of arithmetic precision. No fundamentally new
596 // ideas are involved, but the details are a bit tricky so they are
597 // explained here.
598 //
599 // We have two intertwined needs that lead us to compute bounds on
600 // the prices nodes can have during the assignment computation, on
601 // the assumption that the given problem is feasible:
602 // 1. Infeasibility detection: Infeasibility is detected by
603 // observing that some node's price has been reduced too much by
604 // relabeling operations (see [Goldberg and Tarjan] for the
605 // argument -- duplicated in modified form below -- bounding the
606 // running time of the push/relabel min-cost flow algorithm for
607 // feasible problems); and
608 // 2. Aggressively relabeling nodes and arcs whose matching is
609 // forced: When a left-side node is incident to only one arc a,
610 // any feasible solution must include a, and reducing the price
611 // of Head(a) by any nonnegative amount preserves epsilon-
612 // optimality. Because of this freedom, we'll call this sort of
613 // relabeling (i.e., a relabeling of a right-side node that is
614 // the only neighbor of the left-side node to which it has been
615 // matched in the present double-push operation) a "slack"
616 // relabeling. Relabelings that are not slack relabelings are
617 // called "confined" relabelings. By relabeling Head(a) to have
618 // p(Head(a))=-infinity, we could guarantee that a never becomes
619 // unmatched during the current iteration, and this would prevent
620 // our wasting time repeatedly unmatching and rematching a. But
621 // there are some details we need to handle:
622 // a. The CostValue type cannot represent -infinity;
623 // b. Low node prices are precisely the signal we use to detect
624 // infeasibility (see (1)), so we must be careful not to
625 // falsely conclude that the problem is infeasible as a result
626 // of the low price we gave Head(a); and
627 // c. We need to indicate accurately to the client when our best
628 // understanding indicates that we can't rule out arithmetic
629 // overflow in our calculations. Most importantly, if we don't
630 // warn the client, we must be certain to avoid overflow. This
631 // means our slack relabelings must not be so aggressive as to
632 // create the possibility of unforeseen overflow. Although we
633 // will not achieve this in practice, slack relabelings would
634 // ideally not introduce overflow unless overflow was
635 // inevitable were even the smallest reasonable price change
636 // (== epsilon) used for slack relabelings.
637 // Using the analysis below, we choose a finite amount of price
638 // change for slack relabelings aggressive enough that we don't
639 // waste time doing repeated slack relabelings in a single
640 // iteration, yet modest enough that we keep a good handle on
641 // arithmetic precision and our ability to detect infeasible
642 // problems.
643 //
644 // To provide faithful detection of infeasibility, a dependable
645 // guarantee of arithmetic precision whenever possible, and good
646 // performance by aggressively relabeling nodes whose matching is
647 // forced, we exploit these facts:
648 // 1. Beyond the first iteration, infeasibility detection isn't needed
649 // because a problem is feasible in some iteration if and only if
650 // it's feasible in all others. Therefore we are free to use an
651 // infeasibility detection mechanism that might work in just one
652 // iteration and switch it off in all other iterations.
653 // 2. When we do a slack relabeling, we must choose the amount of
654 // price reduction to use. We choose an amount large enough to
655 // guarantee putting the node's matching to rest, yet (although
656 // we don't bother to prove this explicitly) small enough that
657 // the node's price obeys the overall lower bound that holds if
658 // the slack relabeling amount is small.
659 //
660 // We will establish Propositions (1) and (2) above according to the
661 // following steps:
662 // First, we prove Lemma 1, which is a modified form of lemma 5.8 of
663 // [Goldberg and Tarjan] giving a bound on the difference in price
664 // between the end nodes of certain paths in the residual graph.
665 // Second, we prove Lemma 2, which is technical lemma to establish
666 // reachability of certain "anchor" nodes in the residual graph from
667 // any node where a relabeling takes place.
668 // Third, we apply the first two lemmas to prove Lemma 3 and Lemma
669 // 4, which give two similar bounds that hold whenever the given
670 // problem is feasible: (for feasibility detection) a bound on the
671 // price of any node we relabel during any iteration (and the first
672 // iteration in particular), and (for arithmetic precision) a bound
673 // on the price of any node we relabel during the entire algorithm.
674 //
675 // Finally, we note that if the whole-algorithm price bound can be
676 // represented precisely by the CostValue type, arithmetic overflow
677 // cannot occur (establishing Proposition 1), and assuming no
678 // overflow occurs during the first iteration, any violation of the
679 // first-iteration price bound establishes infeasibility
680 // (Proposition 2).
681 //
682 // The statement of Lemma 1 is perhaps easier to understand when the
683 // reader knows how it will be used. To wit: In this lemma, f' and
684 // e_0 are the flow and error parameter (epsilon) at the beginning
685 // of the current iteration, while f and e_1 are the current
686 // pseudoflow and error parameter when a relabeling of interest
687 // occurs. Without loss of generality, c is the reduced cost
688 // function at the beginning of the current iteration and p is the
689 // change in prices that has taken place in the current iteration.
690 //
691 // Lemma 1 (a variant of lemma 5.8 from [Goldberg and Tarjan]): Let
692 // f be a pseudoflow and let f' be a flow. Suppose P is a simple
693 // path from right-side node v to right-side node w such that P is
694 // residual with respect to f and reverse(P) is residual with
695 // respect to f'. Further, suppose c is an arc cost function with
696 // respect to which f' is e_0-optimal with the zero price function
697 // and p is a price function with respect to which f is e_1-optimal
698 // with respect to p. Then
699 // p(v) - p(w) >= -(e_0 + e_1) * (n-2)/2. (***)
700 //
701 // Proof: We have c_p(P) = p(v) + c(P) - p(w) and hence
702 // p(v) - p(w) = c_p(P) - c(P).
703 // So we seek a bound on c_p(P) - c(P).
704 // p(v) = c_p(P) - c(P).
705 // Let arc a lie on P, which implies that a is residual with respect
706 // to f and reverse(a) is residual with respect to f'.
707 // Case 1: a is a forward arc. Then by e_1-optimality of f with
708 // respect to p, c_p(a) >= 0 and reverse(a) is residual with
709 // respect to f'. By e_0-optimality of f', c(a) <= e_0. So
710 // c_p(a) - c(a) >= -e_0.
711 // Case 2: a is a reverse arc. Then by e_1-optimality of f with
712 // respect to p, c_p(a) >= -e_1 and reverse(a) is residual
713 // with respect to f'. By e_0-optimality of f', c(a) <= 0.
714 // So
715 // c_p(a) - c(a) >= -e_1.
716 // We assumed v and w are both right-side nodes, so there are at
717 // most n - 2 arcs on the path P, of which at most (n-2)/2 are
718 // forward arcs and at most (n-2)/2 are reverse arcs, so
719 // p(v) - p(w) = c_p(P) - c(P)
720 // >= -(e_0 + e_1) * (n-2)/2. (***)
721 //
722 // Some of the rest of our argument is given as a sketch, omitting
723 // several details. Also elided here are some minor technical issues
724 // related to the first iteration, inasmuch as our arguments assume
725 // on the surface a "previous iteration" that doesn't exist in that
726 // case. The issues are not substantial, just a bit messy.
727 //
728 // Lemma 2 is analogous to lemma 5.7 of [Goldberg and Tarjan], where
729 // they have only relabelings that take place at nodes with excess
730 // while we have only relabelings that take place as part of the
731 // double-push operation at nodes without excess.
732 //
733 // Lemma 2: If the problem is feasible, for any node v with excess,
734 // there exists a path P from v to a node w with deficit such that P
735 // is residual with respect to the current pseudoflow, and
736 // reverse(P) is residual with respect to the flow at the beginning
737 // of the current iteration. (Note that such a path exactly
738 // satisfies the conditions of Lemma 1.)
739 //
740 // Let the bound from Lemma 1 with p(w) = 0 be called B(e_0, e_1),
741 // and let us say that when a slack relabeling of a node v occurs,
742 // we will change the price of v by B(e_0, e_1) such that v tightly
743 // satisfies the bound of Lemma 1. Explicitly, we define
744 // B(e_0, e_1) = -(e_0 + e_1) * (n-2)/2.
745 //
746 // Lemma 1 and Lemma 2 combine to bound the price change during an
747 // iteration for any node with excess. Viewed a different way, Lemma
748 // 1 and Lemma 2 tell us that if epsilon-optimality can be preserved
749 // by changing the price of a node by B(e_0, e_1), that node will
750 // never have excess again during the current iteration unless the
751 // problem is infeasible. This insight gives us an approach to
752 // detect infeasibility (by observing prices on nodes with excess
753 // that violate this bound) and to relabel nodes aggressively enough
754 // to avoid unnecessary future work while we also avoid falsely
755 // concluding the problem is infeasible.
756 //
757 // From Lemma 1 and Lemma 2, and taking into account our knowledge
758 // of the slack relabeling amount, we have Lemma 3.
759 //
760 // Lemma 3: During any iteration, if the given problem is feasible
761 // the price of any node is reduced by less than
762 // -2 * B(e_0, e_1) = (e_0 + e_1) * (n-2).
763 //
764 // Proof: Straightforward, omitted for expedience.
765 //
766 // In the case where e_0 = e_1 * alpha, we can express the bound
767 // just in terms of e_1, the current iteration's value of epsilon_:
768 // B(e_1) = B(e_1 * alpha, e_1) = -(1 + alpha) * e_1 * (n-2)/2,
769 // so we have that p(v) is reduced by less than 2 * B(e_1).
770 //
771 // Because we use truncating division to compute each iteration's error
772 // parameter from that of the previous iteration, it isn't exactly
773 // the case that e_0 = e_1 * alpha as we just assumed. To patch this
774 // up, we can use the observation that
775 // e_1 = floor(e_0 / alpha),
776 // which implies
777 // -e_0 > -(e_1 + 1) * alpha
778 // to rewrite from (***):
779 // p(v) > 2 * B(e_0, e_1) > 2 * B((e_1 + 1) * alpha, e_1)
780 // = 2 * -((e_1 + 1) * alpha + e_1) * (n-2)/2
781 // = 2 * -(1 + alpha) * e_1 * (n-2)/2 - alpha * (n-2)
782 // = 2 * B(e_1) - alpha * (n-2)
783 // = -((1 + alpha) * e_1 + alpha) * (n-2).
784 //
785 // We sum up the bounds for all the iterations to get Lemma 4:
786 //
787 // Lemma 4: If the given problem is feasible, after k iterations the
788 // price of any node is always greater than
789 // -((1 + alpha) * C + (k * alpha)) * (n-2)
790 //
791 // Proof: Suppose the price decrease of every node in the iteration
792 // with epsilon_ == x is bounded by B(x) which is proportional to x
793 // (not surprisingly, this will be the same function B() as
794 // above). Assume for simplicity that C, the largest cost magnitude,
795 // is a power of alpha. Then the price of each node, tallied across
796 // all iterations is bounded
797 // p(v) > 2 * B(C/alpha) + 2 * B(C/alpha^2) + ... + 2 * B(kMinEpsilon)
798 // == 2 * B(C/alpha) * alpha / (alpha - 1)
799 // == 2 * B(C) / (alpha - 1).
800 // As above, this needs some patching up to handle the fact that we
801 // use truncating arithmetic. We saw that each iteration effectively
802 // reduces the price bound by alpha * (n-2), hence if there are k
803 // iterations, the bound is
804 // p(v) > 2 * B(C) / (alpha - 1) - k * alpha * (n-2)
805 // = -(1 + alpha) * C * (n-2) / (alpha - 1) - k * alpha * (n-2)
806 // = (n-2) * (C * (1 + alpha) / (1 - alpha) - k * alpha).
807 //
808 // The bound of lemma 4 can be used to warn for possible overflow of
809 // arithmetic precision. But because it involves the number of
810 // iterations, k, we might as well count through the iterations
811 // simply adding up the bounds given by Lemma 3 to get a tighter
812 // result. This is what the implementation does.
813
814 // A lower bound on the price of any node at any time throughout the
815 // computation. A price below this level proves infeasibility; this
816 // value is used for feasibility detection. We use this value also
817 // to rule out the possibility of arithmetic overflow or warn the
818 // client that we have not been able to rule out that possibility.
819 //
820 // We can use the value implied by Lemma 4 here, but note that that
821 // value includes k, the number of iterations. It's plenty fast if
822 // we count through the iterations to compute that value, but if
823 // we're going to count through the iterations, we might as well use
824 // the two-parameter bound from Lemma 3, summing up as we go. This
825 // gives us a tighter bound and more comprehensible code.
826 //
827 // While computing this bound, if we find the value justified by the
828 // theory lies outside the representable range of CostValue, we
829 // conclude that the given arc costs have magnitudes so large that
830 // we cannot guarantee our calculations don't overflow. If the value
831 // justified by the theory lies inside the representable range of
832 // CostValue, we commit that our calculation will not overflow. This
833 // commitment means we need to be careful with the amount by which
834 // we relabel right-side nodes that are incident to any node with
835 // only one neighbor.
836 CostValue price_lower_bound_;
837
838 // A bound on the amount by which a node's price can be reduced
839 // during the current iteration, used only for slack
840 // relabelings. Where epsilon is the first iteration's error
841 // parameter and C is the largest magnitude of an arc cost, we set
842 // slack_relabeling_price_ = -B(C, epsilon)
843 // = (C + epsilon) * (n-2)/2.
844 //
845 // We could use slack_relabeling_price_ for feasibility detection
846 // but the feasibility threshold is double the slack relabeling
847 // amount and we judge it not to be worth having to multiply by two
848 // gratuitously to check feasibility in each double push
849 // operation. Instead we settle for feasibility detection using
850 // price_lower_bound_ instead, which is somewhat slower in the
851 // infeasible case because more relabelings will be required for
852 // some node price to attain the looser bound.
853 CostValue slack_relabeling_price_;
854
855 // Computes the value of the bound on price reduction for an
856 // iteration, given the old and new values of epsilon_. Because the
857 // expression computed here is used in at least one place where we
858 // want an additional factor in the denominator, we take that factor
859 // as an argument. If extra_divisor == 1, this function computes of
860 // the function B() discussed above.
861 //
862 // Avoids overflow in computing the bound, and sets *in_range =
863 // false if the value of the bound doesn't fit in CostValue.
864 inline CostValue PriceChangeBound(CostValue old_epsilon,
865 CostValue new_epsilon,
866 bool* in_range) const {
867 const CostValue n = graph_->num_nodes();
868 // We work in double-precision floating point to determine whether
869 // we'll overflow the integral CostValue type's range of
870 // representation. Switching between integer and double is a
871 // rather expensive operation, but we do this only twice per
872 // scaling iteration, so we can afford it rather than resort to
873 // complex and subtle tricks within the bounds of integer
874 // arithmetic.
875 //
876 // You will want to read the comments above about
877 // price_lower_bound_ and slack_relabeling_price_, and have a
878 // pencil handy. :-)
879 const double result =
880 static_cast<double>(std::max<CostValue>(1, n / 2 - 1)) *
881 (static_cast<double>(old_epsilon) + static_cast<double>(new_epsilon));
882 const double limit =
883 static_cast<double>(std::numeric_limits<CostValue>::max());
884 if (result > limit) {
885 // Our integer computations could overflow.
886 if (in_range != nullptr) *in_range = false;
888 } else {
889 // Don't touch *in_range; other computations could already have
890 // set it to false and we don't want to overwrite that result.
891 return static_cast<CostValue>(result);
892 }
893 }
894
895 // A scaled record of the largest arc-cost magnitude we've been
896 // given during problem setup. This is used to set the initial value
897 // of epsilon_, which in turn is used not only as the error
898 // parameter but also to determine whether we risk arithmetic
899 // overflow during the algorithm.
900 //
901 // Note: Our treatment of arithmetic overflow assumes the following
902 // property of CostValue:
903 // -std::numeric_limits<CostValue>::max() is a representable
904 // CostValue.
905 // That property is satisfied if CostValue uses a two's-complement
906 // representation.
907 CostValue largest_scaled_cost_magnitude_;
908
909 // The total excess in the graph. Given our asymmetric definition of
910 // epsilon-optimality and our use of the double-push operation, this
911 // equals the number of unmatched left-side nodes.
912 NodeIndex total_excess_;
913
914 // Indexed by node index, the price_ values are maintained only for
915 // right-side nodes.
916 //
917 // Note: We use a ZVector to only allocate a vector of size num_left_nodes_
918 // instead of 2*num_left_nodes_ since the right-side node indices start at
919 // num_left_nodes_.
920 ZVector<CostValue> price_;
921
922 // Indexed by left-side node index, the matched_arc_ array gives the
923 // arc index of the arc matching any given left-side node, or
924 // GraphType::kNilArc if the node is unmatched.
925 std::vector<ArcIndex> matched_arc_;
926
927 // Indexed by right-side node index, the matched_node_ array gives
928 // the node index of the left-side node matching any given
929 // right-side node, or GraphType::kNilNode if the right-side node is
930 // unmatched.
931 //
932 // Note: We use a ZVector for the same reason as for price_.
933 ZVector<NodeIndex> matched_node_;
934
935 // The array of arc costs as given in the problem definition, except
936 // that they are scaled up by the number of nodes in the graph so we
937 // can use integer arithmetic throughout.
938 std::vector<CostValue> scaled_arc_cost_;
939
940 // The container of active nodes (i.e., unmatched nodes). This can
941 // be switched easily between ActiveNodeStack and ActiveNodeQueue
942 // for experimentation.
943 std::unique_ptr<ActiveNodeContainerInterface> active_nodes_;
944
945 // Statistics giving the overall numbers of various operations the
946 // algorithm performs.
947 Stats total_stats_;
948
949 // Statistics giving the numbers of various operations the algorithm
950 // has performed in the current iteration.
951 Stats iteration_stats_;
952
953 DISALLOW_COPY_AND_ASSIGN(LinearSumAssignment);
954};
955
956// Implementation of out-of-line LinearSumAssignment template member
957// functions.
958
959template <typename GraphType>
960const CostValue LinearSumAssignment<GraphType>::kMinEpsilon = 1;
961
962template <typename GraphType>
964 const GraphType& graph, const NodeIndex num_left_nodes)
965 : graph_(&graph),
966 num_left_nodes_(num_left_nodes),
967 success_(false),
968 cost_scaling_factor_(1 + num_left_nodes),
969 alpha_(absl::GetFlag(FLAGS_assignment_alpha)),
970 epsilon_(0),
971 price_lower_bound_(0),
972 slack_relabeling_price_(0),
973 largest_scaled_cost_magnitude_(0),
974 total_excess_(0),
975 price_(num_left_nodes, 2 * num_left_nodes - 1),
976 matched_arc_(num_left_nodes, 0),
977 matched_node_(num_left_nodes, 2 * num_left_nodes - 1),
978 scaled_arc_cost_(graph.max_end_arc_index(), 0),
979 active_nodes_(absl::GetFlag(FLAGS_assignment_stack_order)
980 ? static_cast<ActiveNodeContainerInterface*>(
981 new ActiveNodeStack())
982 : static_cast<ActiveNodeContainerInterface*>(
983 new ActiveNodeQueue())) {}
984
985template <typename GraphType>
987 const NodeIndex num_left_nodes, const ArcIndex num_arcs)
988 : graph_(nullptr),
989 num_left_nodes_(num_left_nodes),
990 success_(false),
991 cost_scaling_factor_(1 + num_left_nodes),
992 alpha_(absl::GetFlag(FLAGS_assignment_alpha)),
993 epsilon_(0),
994 price_lower_bound_(0),
995 slack_relabeling_price_(0),
996 largest_scaled_cost_magnitude_(0),
997 total_excess_(0),
998 price_(num_left_nodes, 2 * num_left_nodes - 1),
999 matched_arc_(num_left_nodes, 0),
1000 matched_node_(num_left_nodes, 2 * num_left_nodes - 1),
1001 scaled_arc_cost_(num_arcs, 0),
1002 active_nodes_(absl::GetFlag(FLAGS_assignment_stack_order)
1003 ? static_cast<ActiveNodeContainerInterface*>(
1004 new ActiveNodeStack())
1005 : static_cast<ActiveNodeContainerInterface*>(
1006 new ActiveNodeQueue())) {}
1007
1008template <typename GraphType>
1010 if (graph_ != nullptr) {
1011 DCHECK_GE(arc, 0);
1012 DCHECK_LT(arc, graph_->num_arcs());
1013 NodeIndex head = Head(arc);
1014 DCHECK_LE(num_left_nodes_, head);
1015 }
1016 cost *= cost_scaling_factor_;
1017 const CostValue cost_magnitude = std::abs(cost);
1018 largest_scaled_cost_magnitude_ =
1019 std::max(largest_scaled_cost_magnitude_, cost_magnitude);
1020 scaled_arc_cost_[arc] = cost;
1021}
1022
1023template <typename ArcIndexType>
1025 public:
1026 explicit CostValueCycleHandler(std::vector<CostValue>* cost)
1027 : temp_(0), cost_(cost) {}
1028
1029 void SetTempFromIndex(ArcIndexType source) override {
1030 temp_ = (*cost_)[source];
1031 }
1032
1033 void SetIndexFromIndex(ArcIndexType source,
1034 ArcIndexType destination) const override {
1035 (*cost_)[destination] = (*cost_)[source];
1036 }
1037
1038 void SetIndexFromTemp(ArcIndexType destination) const override {
1039 (*cost_)[destination] = temp_;
1040 }
1041
1043
1044 private:
1045 CostValue temp_;
1046 std::vector<CostValue>* const cost_;
1047
1048 DISALLOW_COPY_AND_ASSIGN(CostValueCycleHandler);
1049};
1050
1051// Logically this class should be defined inside OptimizeGraphLayout,
1052// but compilation fails if we do that because C++98 doesn't allow
1053// instantiation of member templates with function-scoped types as
1054// template parameters, which in turn is because those function-scoped
1055// types lack linkage.
1056template <typename GraphType>
1058 public:
1059 explicit ArcIndexOrderingByTailNode(const GraphType& graph) : graph_(graph) {}
1060
1061 // Says ArcIndex a is less than ArcIndex b if arc a's tail is less
1062 // than arc b's tail. If their tails are equal, orders according to
1063 // heads.
1065 typename GraphType::ArcIndex b) const {
1066 return ((graph_.Tail(a) < graph_.Tail(b)) ||
1067 ((graph_.Tail(a) == graph_.Tail(b)) &&
1068 (graph_.Head(a) < graph_.Head(b))));
1069 }
1070
1071 private:
1072 const GraphType& graph_;
1073
1074 // Copy and assign are allowed; they have to be for STL to work
1075 // with this functor, although it seems like a bug for STL to be
1076 // written that way.
1077};
1078
1079// Passes ownership of the cycle handler to the caller.
1080template <typename GraphType>
1081PermutationCycleHandler<typename GraphType::ArcIndex>*
1084 &scaled_arc_cost_);
1085}
1086
1087template <typename GraphType>
1089 // The graph argument is only to give us a non-const-qualified
1090 // handle on the graph we already have. Any different graph is
1091 // nonsense.
1092 DCHECK_EQ(graph_, graph);
1093 const ArcIndexOrderingByTailNode<GraphType> compare(*graph_);
1095 &scaled_arc_cost_);
1096 TailArrayManager<GraphType> tail_array_manager(graph);
1098 graph->GroupForwardArcsByFunctor(compare, &cycle_handler);
1099 tail_array_manager.ReleaseTailArrayIfForwardGraph();
1100}
1101
1102template <typename GraphType>
1104 const CostValue current_epsilon) const {
1105 return std::max(current_epsilon / alpha_, kMinEpsilon);
1106}
1107
1108template <typename GraphType>
1109bool LinearSumAssignment<GraphType>::UpdateEpsilon() {
1110 CostValue new_epsilon = NewEpsilon(epsilon_);
1111 slack_relabeling_price_ = PriceChangeBound(epsilon_, new_epsilon, nullptr);
1112 epsilon_ = new_epsilon;
1113 VLOG(3) << "Updated: epsilon_ == " << epsilon_;
1114 VLOG(4) << "slack_relabeling_price_ == " << slack_relabeling_price_;
1115 DCHECK_GT(slack_relabeling_price_, 0);
1116 // For today we always return true; in the future updating epsilon
1117 // in sophisticated ways could conceivably detect infeasibility
1118 // before the first iteration of Refine().
1119 return true;
1120}
1121
1122// For production code that checks whether a left-side node is active.
1123template <typename GraphType>
1124inline bool LinearSumAssignment<GraphType>::IsActive(
1125 NodeIndex left_node) const {
1126 DCHECK_LT(left_node, num_left_nodes_);
1127 return matched_arc_[left_node] == GraphType::kNilArc;
1128}
1129
1130// Only for debugging. Separate from the production IsActive() method
1131// so that method can assert that its argument is a left-side node,
1132// while for debugging we need to be able to test any node.
1133template <typename GraphType>
1134inline bool LinearSumAssignment<GraphType>::IsActiveForDebugging(
1135 NodeIndex node) const {
1136 if (node < num_left_nodes_) {
1137 return IsActive(node);
1138 } else {
1139 return matched_node_[node] == GraphType::kNilNode;
1140 }
1141}
1142
1143template <typename GraphType>
1144void LinearSumAssignment<GraphType>::InitializeActiveNodeContainer() {
1145 DCHECK(active_nodes_->Empty());
1146 for (BipartiteLeftNodeIterator node_it(*graph_, num_left_nodes_);
1147 node_it.Ok(); node_it.Next()) {
1148 const NodeIndex node = node_it.Index();
1149 if (IsActive(node)) {
1150 active_nodes_->Add(node);
1151 }
1152 }
1153}
1154
1155// There exists a price function such that the admissible arcs at the
1156// beginning of an iteration are exactly the reverse arcs of all
1157// matching arcs. Saturating all admissible arcs with respect to that
1158// price function therefore means simply unmatching every matched
1159// node.
1160//
1161// In the future we will price out arcs, which will reduce the set of
1162// nodes we unmatch here. If a matching arc is priced out, we will not
1163// unmatch its endpoints since that element of the matching is
1164// guaranteed not to change.
1165template <typename GraphType>
1166void LinearSumAssignment<GraphType>::SaturateNegativeArcs() {
1167 total_excess_ = 0;
1168 for (BipartiteLeftNodeIterator node_it(*graph_, num_left_nodes_);
1169 node_it.Ok(); node_it.Next()) {
1170 const NodeIndex node = node_it.Index();
1171 if (IsActive(node)) {
1172 // This can happen in the first iteration when nothing is
1173 // matched yet.
1174 total_excess_ += 1;
1175 } else {
1176 // We're about to create a unit of excess by unmatching these nodes.
1177 total_excess_ += 1;
1178 const NodeIndex mate = GetMate(node);
1179 matched_arc_[node] = GraphType::kNilArc;
1180 matched_node_[mate] = GraphType::kNilNode;
1181 }
1182 }
1183}
1184
1185// Returns true for success, false for infeasible.
1186template <typename GraphType>
1187bool LinearSumAssignment<GraphType>::DoublePush(NodeIndex source) {
1188 DCHECK_GT(num_left_nodes_, source);
1189 DCHECK(IsActive(source)) << "Node " << source
1190 << "must be active (unmatched)!";
1191 ImplicitPriceSummary summary = BestArcAndGap(source);
1192 const ArcIndex best_arc = summary.first;
1193 const CostValue gap = summary.second;
1194 // Now we have the best arc incident to source, i.e., the one with
1195 // minimum reduced cost. Match that arc, unmatching its head if
1196 // necessary.
1197 if (best_arc == GraphType::kNilArc) {
1198 return false;
1199 }
1200 const NodeIndex new_mate = Head(best_arc);
1201 const NodeIndex to_unmatch = matched_node_[new_mate];
1202 if (to_unmatch != GraphType::kNilNode) {
1203 // Unmatch new_mate from its current mate, pushing the unit of
1204 // flow back to a node on the left side as a unit of excess.
1205 matched_arc_[to_unmatch] = GraphType::kNilArc;
1206 active_nodes_->Add(to_unmatch);
1207 // This counts as a double push.
1208 iteration_stats_.double_pushes_ += 1;
1209 } else {
1210 // We are about to increase the cardinality of the matching.
1211 total_excess_ -= 1;
1212 // This counts as a single push.
1213 iteration_stats_.pushes_ += 1;
1214 }
1215 matched_arc_[source] = best_arc;
1216 matched_node_[new_mate] = source;
1217 // Finally, relabel new_mate.
1218 iteration_stats_.relabelings_ += 1;
1219 const CostValue new_price = price_[new_mate] - gap - epsilon_;
1220 price_[new_mate] = new_price;
1221 return new_price >= price_lower_bound_;
1222}
1223
1224template <typename GraphType>
1225bool LinearSumAssignment<GraphType>::Refine() {
1226 SaturateNegativeArcs();
1227 InitializeActiveNodeContainer();
1228 while (total_excess_ > 0) {
1229 // Get an active node (i.e., one with excess == 1) and discharge
1230 // it using DoublePush.
1231 const NodeIndex node = active_nodes_->Get();
1232 if (!DoublePush(node)) {
1233 // Infeasibility detected.
1234 //
1235 // If infeasibility is detected after the first iteration, we
1236 // have a bug. We don't crash production code in this case but
1237 // we know we're returning a wrong answer so we we leave a
1238 // message in the logs to increase our hope of chasing down the
1239 // problem.
1240 LOG_IF(DFATAL, total_stats_.refinements_ > 0)
1241 << "Infeasibility detection triggered after first iteration found "
1242 << "a feasible assignment!";
1243 return false;
1244 }
1245 }
1246 DCHECK(active_nodes_->Empty());
1247 iteration_stats_.refinements_ += 1;
1248 return true;
1249}
1250
1251// Computes best_arc, the minimum reduced-cost arc incident to
1252// left_node and admissibility_gap, the amount by which the reduced
1253// cost of best_arc must be increased to make it equal in reduced cost
1254// to another residual arc incident to left_node.
1255//
1256// Precondition: left_node is unmatched and has at least one incident
1257// arc. This allows us to simplify the code. The debug-only
1258// counterpart to this routine is LinearSumAssignment::ImplicitPrice()
1259// and it assumes there is an incident arc but does not assume the
1260// node is unmatched. The condition that each left node has at least
1261// one incident arc is explicitly computed during FinalizeSetup().
1262//
1263// This function is large enough that our suggestion that the compiler
1264// inline it might be pointless.
1265template <typename GraphType>
1266inline typename LinearSumAssignment<GraphType>::ImplicitPriceSummary
1267LinearSumAssignment<GraphType>::BestArcAndGap(NodeIndex left_node) const {
1268 DCHECK(IsActive(left_node))
1269 << "Node " << left_node << " must be active (unmatched)!";
1270 DCHECK_GT(epsilon_, 0);
1271 typename GraphType::OutgoingArcIterator arc_it(*graph_, left_node);
1272 ArcIndex best_arc = arc_it.Index();
1273 CostValue min_partial_reduced_cost = PartialReducedCost(best_arc);
1274 // We choose second_min_partial_reduced_cost so that in the case of
1275 // the largest possible gap (which results from a left-side node
1276 // with only a single incident residual arc), the corresponding
1277 // right-side node will be relabeled by an amount that exactly
1278 // matches slack_relabeling_price_.
1279 const CostValue max_gap = slack_relabeling_price_ - epsilon_;
1280 CostValue second_min_partial_reduced_cost =
1281 min_partial_reduced_cost + max_gap;
1282 for (arc_it.Next(); arc_it.Ok(); arc_it.Next()) {
1283 const ArcIndex arc = arc_it.Index();
1284 const CostValue partial_reduced_cost = PartialReducedCost(arc);
1285 if (partial_reduced_cost < second_min_partial_reduced_cost) {
1286 if (partial_reduced_cost < min_partial_reduced_cost) {
1287 best_arc = arc;
1288 second_min_partial_reduced_cost = min_partial_reduced_cost;
1289 min_partial_reduced_cost = partial_reduced_cost;
1290 } else {
1291 second_min_partial_reduced_cost = partial_reduced_cost;
1292 }
1293 }
1294 }
1295 const CostValue gap = std::min<CostValue>(
1296 second_min_partial_reduced_cost - min_partial_reduced_cost, max_gap);
1297 DCHECK_GE(gap, 0);
1298 return std::make_pair(best_arc, gap);
1299}
1300
1301// Only for debugging.
1302//
1303// Requires the precondition, explicitly computed in FinalizeSetup(),
1304// that every left-side node has at least one incident arc.
1305template <typename GraphType>
1306inline CostValue LinearSumAssignment<GraphType>::ImplicitPrice(
1307 NodeIndex left_node) const {
1308 DCHECK_GT(num_left_nodes_, left_node);
1309 DCHECK_GT(epsilon_, 0);
1310 typename GraphType::OutgoingArcIterator arc_it(*graph_, left_node);
1311 // We must not execute this method if left_node has no incident arc.
1312 DCHECK(arc_it.Ok());
1313 ArcIndex best_arc = arc_it.Index();
1314 if (best_arc == matched_arc_[left_node]) {
1315 arc_it.Next();
1316 if (arc_it.Ok()) {
1317 best_arc = arc_it.Index();
1318 }
1319 }
1320 CostValue min_partial_reduced_cost = PartialReducedCost(best_arc);
1321 if (!arc_it.Ok()) {
1322 // Only one arc is incident to left_node, and the node is
1323 // currently matched along that arc, which must be the case in any
1324 // feasible solution. Therefore we implicitly price this node so
1325 // low that we will never consider unmatching it.
1326 return -(min_partial_reduced_cost + slack_relabeling_price_);
1327 }
1328 for (arc_it.Next(); arc_it.Ok(); arc_it.Next()) {
1329 const ArcIndex arc = arc_it.Index();
1330 if (arc != matched_arc_[left_node]) {
1331 const CostValue partial_reduced_cost = PartialReducedCost(arc);
1332 if (partial_reduced_cost < min_partial_reduced_cost) {
1333 min_partial_reduced_cost = partial_reduced_cost;
1334 }
1335 }
1336 }
1337 return -min_partial_reduced_cost;
1338}
1339
1340// Only for debugging.
1341template <typename GraphType>
1342bool LinearSumAssignment<GraphType>::AllMatched() const {
1343 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
1344 if (IsActiveForDebugging(node)) {
1345 return false;
1346 }
1347 }
1348 return true;
1349}
1350
1351// Only for debugging.
1352template <typename GraphType>
1353bool LinearSumAssignment<GraphType>::EpsilonOptimal() const {
1354 for (BipartiteLeftNodeIterator node_it(*graph_, num_left_nodes_);
1355 node_it.Ok(); node_it.Next()) {
1356 const NodeIndex left_node = node_it.Index();
1357 // Get the implicit price of left_node and make sure the reduced
1358 // costs of left_node's incident arcs are in bounds.
1359 CostValue left_node_price = ImplicitPrice(left_node);
1360 for (typename GraphType::OutgoingArcIterator arc_it(*graph_, left_node);
1361 arc_it.Ok(); arc_it.Next()) {
1362 const ArcIndex arc = arc_it.Index();
1363 const CostValue reduced_cost = left_node_price + PartialReducedCost(arc);
1364 // Note the asymmetric definition of epsilon-optimality that we
1365 // use because it means we can saturate all admissible arcs in
1366 // the beginning of Refine() just by unmatching all matched
1367 // nodes.
1368 if (matched_arc_[left_node] == arc) {
1369 // The reverse arc is residual. Epsilon-optimality requires
1370 // that the reduced cost of the forward arc be at most
1371 // epsilon_.
1372 if (reduced_cost > epsilon_) {
1373 return false;
1374 }
1375 } else {
1376 // The forward arc is residual. Epsilon-optimality requires
1377 // that the reduced cost of the forward arc be at least zero.
1378 if (reduced_cost < 0) {
1379 return false;
1380 }
1381 }
1382 }
1383 }
1384 return true;
1385}
1386
1387template <typename GraphType>
1389 incidence_precondition_satisfied_ = true;
1390 // epsilon_ must be greater than kMinEpsilon so that in the case
1391 // where the largest arc cost is zero, we still do a Refine()
1392 // iteration.
1393 epsilon_ = std::max(largest_scaled_cost_magnitude_, kMinEpsilon + 1);
1394 VLOG(2) << "Largest given cost magnitude: "
1395 << largest_scaled_cost_magnitude_ / cost_scaling_factor_;
1396 // Initialize left-side node-indexed arrays and check incidence
1397 // precondition.
1398 for (NodeIndex node = 0; node < num_left_nodes_; ++node) {
1399 matched_arc_[node] = GraphType::kNilArc;
1400 typename GraphType::OutgoingArcIterator arc_it(*graph_, node);
1401 if (!arc_it.Ok()) {
1402 incidence_precondition_satisfied_ = false;
1403 }
1404 }
1405 // Initialize right-side node-indexed arrays. Example: prices are
1406 // stored only for right-side nodes.
1407 for (NodeIndex node = num_left_nodes_; node < graph_->num_nodes(); ++node) {
1408 price_[node] = 0;
1409 matched_node_[node] = GraphType::kNilNode;
1410 }
1411 bool in_range = true;
1412 double double_price_lower_bound = 0.0;
1413 CostValue new_error_parameter;
1414 CostValue old_error_parameter = epsilon_;
1415 do {
1416 new_error_parameter = NewEpsilon(old_error_parameter);
1417 double_price_lower_bound -=
1418 2.0 * static_cast<double>(PriceChangeBound(
1419 old_error_parameter, new_error_parameter, &in_range));
1420 old_error_parameter = new_error_parameter;
1421 } while (new_error_parameter != kMinEpsilon);
1422 const double limit =
1423 -static_cast<double>(std::numeric_limits<CostValue>::max());
1424 if (double_price_lower_bound < limit) {
1425 in_range = false;
1426 price_lower_bound_ = -std::numeric_limits<CostValue>::max();
1427 } else {
1428 price_lower_bound_ = static_cast<CostValue>(double_price_lower_bound);
1429 }
1430 VLOG(4) << "price_lower_bound_ == " << price_lower_bound_;
1431 DCHECK_LE(price_lower_bound_, 0);
1432 if (!in_range) {
1433 LOG(WARNING) << "Price change bound exceeds range of representable "
1434 << "costs; arithmetic overflow is not ruled out and "
1435 << "infeasibility might go undetected.";
1436 }
1437 return in_range;
1438}
1439
1440template <typename GraphType>
1442 total_stats_.Add(iteration_stats_);
1443 VLOG(3) << "Iteration stats: " << iteration_stats_.StatsString();
1444 iteration_stats_.Clear();
1445}
1446
1447template <typename GraphType>
1449 CHECK(graph_ != nullptr);
1450 bool ok = graph_->num_nodes() == 2 * num_left_nodes_;
1451 if (!ok) return false;
1452 // Note: FinalizeSetup() might have been called already by white-box
1453 // test code or by a client that wants to react to the possibility
1454 // of overflow before solving the given problem, but FinalizeSetup()
1455 // is idempotent and reasonably fast, so we call it unconditionally
1456 // here.
1457 FinalizeSetup();
1458 ok = ok && incidence_precondition_satisfied_;
1459 DCHECK(!ok || EpsilonOptimal());
1460 while (ok && epsilon_ > kMinEpsilon) {
1461 ok = ok && UpdateEpsilon();
1462 ok = ok && Refine();
1463 ReportAndAccumulateStats();
1464 DCHECK(!ok || EpsilonOptimal());
1465 DCHECK(!ok || AllMatched());
1466 }
1467 success_ = ok;
1468 VLOG(1) << "Overall stats: " << total_stats_.StatsString();
1469 return ok;
1470}
1471
1472template <typename GraphType>
1474 // It is illegal to call this method unless we successfully computed
1475 // an optimum assignment.
1476 DCHECK(success_);
1477 CostValue cost = 0;
1478 for (BipartiteLeftNodeIterator node_it(*this); node_it.Ok(); node_it.Next()) {
1479 cost += GetAssignmentCost(node_it.Index());
1480 }
1481 return cost;
1482}
1483
1484} // namespace operations_research
1485
1486#endif // OR_TOOLS_GRAPH_LINEAR_ASSIGNMENT_H_
int64 max
Definition: alldiff_cst.cc:139
#define LOG_IF(severity, condition)
Definition: base/logging.h:479
#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 DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK_GT(val1, val2)
Definition: base/logging.h:890
#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
bool operator()(typename GraphType::ArcIndex a, typename GraphType::ArcIndex b) const
void SetIndexFromIndex(ArcIndexType source, ArcIndexType destination) const override
CostValueCycleHandler(std::vector< CostValue > *cost)
void SetTempFromIndex(ArcIndexType source) override
void SetIndexFromTemp(ArcIndexType destination) const override
BipartiteLeftNodeIterator(const GraphType &graph, NodeIndex num_left_nodes)
ArcIndex GetAssignmentArc(NodeIndex left_node) const
void SetArcCost(ArcIndex arc, CostValue cost)
NodeIndex Head(ArcIndex arc) const
CostValue ArcCost(ArcIndex arc) const
CostValue GetAssignmentCost(NodeIndex node) const
operations_research::PermutationCycleHandler< typename GraphType::ArcIndex > * ArcAnnotationCycleHandler()
NodeIndex GetMate(NodeIndex left_node) const
void SetGraph(const GraphType *graph)
LinearSumAssignment(const GraphType &graph, NodeIndex num_left_nodes)
bool BuildTailArrayFromAdjacencyListsIfForwardGraph() const
Definition: ebert_graph.h:1920
int64_t int64
ABSL_DECLARE_FLAG(int64, assignment_alpha)
const int WARNING
Definition: log_severity.h:31
Definition: cleanup.h:22
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 cost
int64 head