Generated on: Thu Mar 29 07:46:58 PDT 2012 for custom file set | ||
|
||
00001 // Copyright 2010-2012 Google 00002 // Licensed under the Apache License, Version 2.0 (the "License"); 00003 // you may not use this file except in compliance with the License. 00004 // You may obtain a copy of the License at 00005 // 00006 // http://www.apache.org/licenses/LICENSE-2.0 00007 // 00008 // Unless required by applicable law or agreed to in writing, software 00009 // distributed under the License is distributed on an "AS IS" BASIS, 00010 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 00011 // See the License for the specific language governing permissions and 00012 // limitations under the License. 00013 // 00014 // An implementation of a cost-scaling push-relabel algorithm for the 00015 // assignment problem (minimum-cost perfect bipartite matching), from 00016 // the paper of Goldberg and Kennedy (1995). 00017 // 00018 // This implementation finds the minimum-cost perfect assignment in 00019 // the given graph with integral edge weights set through the 00020 // SetArcCost method. 00021 // 00022 // Example usage: 00023 // 00024 // #include "graph/ebert_graph.h" 00025 // #include "graph/linear_assignment.h" 00026 // ... 00027 // ::operations_research::NodeIndex num_nodes = ...; 00028 // ::operations_research::NodeIndex num_left_nodes = num_nodes / 2; 00029 // // Define a num_nodes/2 by num_nodes/2 assignment problem: 00030 // ::operations_research::ArcIndex num_forward_arcs = ...; 00031 // ::operations_research::ForwardStarGraph g(num_nodes, num_arcs); 00032 // ::operations_research::LinearSumAssignment< 00033 // ::operations_research::ForwardStarGraph> a(g, num_left_nodes); 00034 // for (int i = 0; i < num_forward_arcs; ++i) { 00035 // ::operations_research::NodeIndex this_arc_head = ...; 00036 // ::operations_research::NodeIndex this_arc_tail = ...; 00037 // ::operations_research::CostValue this_arc_cost = ...; 00038 // ::operations_research::ArcIndex this_arc_index = 00039 // g.AddArc(this_arc_tail, this_arc_head); 00040 // a.SetArcCost(this_arc_index, this_arc_cost); 00041 // } 00042 // // Compute the optimum assignment. 00043 // bool success = a.ComputeAssignment(); 00044 // // Retrieve the cost of the optimum assignment. 00045 // CostValue optimum_cost = a.GetCost(); 00046 // // Retrieve the node-node correspondence of the optimum assignment and the 00047 // // cost of each node pairing. 00048 // for (::operations_research::LinearSumAssignment::BipartiteLeftNodeIterator 00049 // node_it(a); 00050 // node_it.Ok(); 00051 // node_it.Next()) { 00052 // ::operations_research::NodeIndex left_node = node_it.Index(); 00053 // ::operations_research::NodeIndex right_node = a.GetMate(left_node); 00054 // ::operations_research::CostValue node_pair_cost = 00055 // a.GetAssignmentCost(left_node); 00056 // ... 00057 // } 00058 // 00059 // In the following, we consider a bipartite graph 00060 // G = (V = X union Y, E subset XxY), 00061 // where V denodes the set of nodes (vertices) in the graph, E denotes 00062 // the set of arcs (edges), n = |V| denotes the number of nodes in the 00063 // graph, and m = |E| denotes the number of arcs in the graph. 00064 // 00065 // The set of nodes is divided into two parts, X and Y, and every arc 00066 // must go between a node of X and a node of Y. With each arc is 00067 // associated a cost c(v, w). A matching M is a subset of E with the 00068 // property that no two arcs in M have a head or tail node in common, 00069 // and a perfect matching is a matching that touches every node in the 00070 // graph. The cost of a matching M is the sum of the costs of all the 00071 // arcs in M. 00072 // 00073 // The assignment problem is to find a perfect matching of minimum 00074 // cost in the given bipartite graph. The present algorithm reduces 00075 // the assignment problem to an instance of the minimum-cost flow 00076 // problem and takes advantage of special properties of the resulting 00077 // minimum-cost flow problem to solve it efficiently using a 00078 // push-relabel method. For more information about minimum-cost flow 00079 // see or-tools/graph/min_cost_flow.h 00080 // 00081 // The method used here is the cost-scaling approach for the 00082 // minimum-cost circulation problem as described in [Goldberg and 00083 // Tarjan] with some technical modifications: 00084 // 1. For efficiency, we solve a transportation problem instead of 00085 // minimum-cost circulation. We might revisit this decision if it 00086 // is important to handle problems in which no perfect matching 00087 // exists. 00088 // 2. We use a modified "asymmetric" notion of epsilon-optimality in 00089 // which left-to-right residual arcs are required to have reduced 00090 // cost bounded below by zero and right-to-left residual arcs are 00091 // required to have reduced cost bounded below by -epsilon. For 00092 // each residual arc direction, the reduced-cost threshold for 00093 // admissibility is epsilon/2 above the threshold for epsilon 00094 // optimality. 00095 // 3. We do not limit the applicability of the relabeling operation to 00096 // nodes with excess. Instead we use the double-push operation 00097 // (discussed in the Goldberg and Kennedy CSA paper and Kennedy's 00098 // thesis) which relabels right-side nodes just *after* they have 00099 // been discharged. 00100 // The above differences are explained in detail in [Kennedy's thesis] 00101 // and explained not quite as cleanly in [Goldberg and Kennedy's CSA 00102 // paper]. But note that the thesis explanation uses a value of 00103 // epsilon that's double what we use here. 00104 // 00105 // Some definitions: 00106 // Active: A node is called active when it has excess. It is 00107 // eligible to be pushed from. In this implementation, every active 00108 // node is on the left side of the graph where prices are determined 00109 // implicitly, so no left-side relabeling is necessary before 00110 // pushing from an active node. We do, however, need to compute 00111 // the implications for price changes on the affected right-side 00112 // nodes. 00113 // Admissible: A residual arc (one that can carry more flow) is 00114 // called admissible when its reduced cost is small enough. We can 00115 // push additional flow along such an arc without violating 00116 // epsilon-optimality. In the case of a left-to-right residual 00117 // arc, the reduced cost must be at most epsilon/2. In the case of 00118 // a right-to-left residual arc, the reduced cost must be at most 00119 // -epsilon/2. The careful reader will note that these thresholds 00120 // are not used explicitly anywhere in this implementation, and 00121 // the reason is the implicit pricing of left-side nodes. 00122 // Reduced cost: Essentially an arc's reduced cost is its 00123 // complementary slackness. In push-relabel algorithms this is 00124 // c_p(v, w) = p(v) + c(v, w) - p(w), 00125 // where p() is the node price function and c(v, w) is the cost of 00126 // the arc from v to w. See min_cost_flow.h for more details. 00127 // Partial reduced cost: We maintain prices implicitly for left-side 00128 // nodes in this implementation, so instead of reduced costs we 00129 // work with partial reduced costs, defined as 00130 // c'_p(v, w) = c(v, w) - p(w). 00131 // 00132 // We check at initialization time for the possibility of arithmetic 00133 // overflow and warn if the given costs are too large. In many cases 00134 // the bound we use to trigger the warning is pessimistic so the given 00135 // problem can often be solved even if we warn that overflow is 00136 // possible. 00137 // 00138 // We don't use the interface from 00139 // operations_research/algorithms/hungarian.h because we want to be 00140 // able to express sparse problems efficiently. 00141 // 00142 // When asked to solve the given assignment problem we return a 00143 // boolean to indicate whether the given problem was feasible. 00144 // 00145 // References: 00146 // [ Goldberg and Kennedy's CSA paper ] A. V. Goldberg and R. Kennedy, 00147 // "An Efficient Cost Scaling Algorithm for the Assignment Problem." 00148 // Mathematical Programming, Vol. 71, pages 153-178, December 1995. 00149 // 00150 // [ Goldberg and Tarjan ] A. V. Goldberg and R. E. Tarjan, "Finding 00151 // Minimum-Cost Circulations by Successive Approximation." Mathematics 00152 // of Operations Research, Vol. 15, No. 3, pages 430-466, August 1990. 00153 // 00154 // [ Kennedy's thesis ] J. R. Kennedy, Jr., "Solving Unweighted and 00155 // Weighted Bipartite Matching Problems in Theory and Practice." 00156 // Stanford University Doctoral Dissertation, Department of Computer 00157 // Science, 1995. 00158 // 00159 // [ Burkard et al. ] R. Burkard, M. Dell'Amico, S. Martello, "Assignment 00160 // Problems", SIAM, 2009, ISBN: 978-0898716634, 00161 // http://www.amazon.com/dp/0898716632/ 00162 // 00163 // [ Ahuja et al. ] R. K. Ahuja, T. L. Magnanti, J. B. Orlin, "Network Flows: 00164 // Theory, Algorithms, and Applications," Prentice Hall, 1993, 00165 // ISBN: 978-0136175490, http://www.amazon.com/dp/013617549X 00166 // 00167 // Keywords: linear sum assignment problem, Hungarian method, Goldberg, Kennedy. 00168 00169 #ifndef OR_TOOLS_GRAPH_LINEAR_ASSIGNMENT_H_ 00170 #define OR_TOOLS_GRAPH_LINEAR_ASSIGNMENT_H_ 00171 00172 #include <algorithm> 00173 #include <cstdlib> 00174 #include <deque> 00175 #include <limits> 00176 #include <string> 00177 #include <utility> 00178 #include <vector> 00179 00180 #include "base/commandlineflags.h" 00181 #include "base/integral_types.h" 00182 #include "base/logging.h" 00183 #include "base/macros.h" 00184 #include "base/scoped_ptr.h" 00185 #include "base/stringprintf.h" 00186 #include "graph/ebert_graph.h" 00187 #include "util/permutation.h" 00188 00189 using std::string; 00190 00191 #ifndef SWIG 00192 DECLARE_int64(assignment_alpha); 00193 DECLARE_int32(assignment_progress_logging_period); 00194 DECLARE_bool(assignment_stack_order); 00195 #endif 00196 00197 namespace operations_research { 00198 00199 template <typename GraphType> class LinearSumAssignment { 00200 public: 00201 #ifndef SWIG 00202 #endif 00203 00204 // This class modifies the given graph by adding arcs to it as costs 00205 // are specified via SetArcCost, but does not take ownership. 00206 LinearSumAssignment(const GraphType& graph, NodeIndex num_left_nodes); 00207 virtual ~LinearSumAssignment() {} 00208 00209 // Sets the cost-scaling divisor, i.e., the amount by which we 00210 // divide the scaling parameter on each iteration. 00211 void SetCostScalingDivisor(CostValue factor) { 00212 alpha_ = factor; 00213 } 00214 00215 // Optimizes the layout of the graph for the access pattern our 00216 // implementation will use. 00217 void OptimizeGraphLayout(GraphType* graph); 00218 00219 // Allows tests, iterators, etc., to inspect our underlying graph. 00220 inline const GraphType& Graph() const { return graph_; } 00221 00222 // These handy member functions make the code more compact, and we 00223 // expose them to clients so that client code that doesn't have 00224 // direct access to the graph can learn about the optimum assignment 00225 // once it is computed. 00226 inline NodeIndex Head(ArcIndex arc) const { 00227 return graph_.Head(arc); 00228 } 00229 00230 // Returns the original arc cost for use by a client that's 00231 // iterating over the optimum assignment. 00232 virtual CostValue ArcCost(ArcIndex arc) const { 00233 DCHECK_EQ(0, scaled_arc_cost_[arc] % cost_scaling_factor_); 00234 return scaled_arc_cost_[arc] / cost_scaling_factor_; 00235 } 00236 00237 // Sets the cost of an arc already present in the given graph. 00238 virtual void SetArcCost(ArcIndex arc, 00239 CostValue cost); 00240 00241 // Completes initialization after the problem is fully specified. 00242 // Returns true if we successfully prove that arithmetic 00243 // calculations are guaranteed not to overflow. ComputeAssignment() 00244 // calls this method itself, so only clients that care about 00245 // obtaining a warning about the possibility of arithmetic precision 00246 // problems need to call this method explicitly. 00247 // 00248 // Separate from ComputeAssignment() for white-box testing and for 00249 // clients that need to react to the possibility that arithmetic 00250 // overflow is not ruled out. 00251 // 00252 // FinalizeSetup() is idempotent. 00253 virtual bool FinalizeSetup(); 00254 00255 // Computes the optimum assignment. Returns true on success. Return 00256 // value of false implies the given problem is infeasible. 00257 virtual bool ComputeAssignment(); 00258 00259 // Returns the cost of the minimum-cost perfect matching. 00260 // Precondition: success_ == true, signifying that we computed the 00261 // optimum assignment for a feasible problem. 00262 virtual CostValue GetCost() const; 00263 00264 // Returns the total number of nodes in the given problem. 00265 virtual NodeIndex NumNodes() const { 00266 return graph_.num_nodes(); 00267 } 00268 00269 // Returns the number of nodes on the left side of the given 00270 // problem. 00271 virtual NodeIndex NumLeftNodes() const { 00272 return num_left_nodes_; 00273 } 00274 00275 // Returns the arc through which the given node is matched. 00276 inline ArcIndex GetAssignmentArc(NodeIndex left_node) const { 00277 DCHECK_LT(left_node, num_left_nodes_); 00278 return matched_arc_[left_node]; 00279 } 00280 00281 // Returns the cost of the assignment arc incident to the given 00282 // node. 00283 inline CostValue GetAssignmentCost(NodeIndex node) const { 00284 return ArcCost(GetAssignmentArc(node)); 00285 } 00286 00287 // Returns the node to which the given node is matched. 00288 inline NodeIndex GetMate(NodeIndex left_node) const { 00289 DCHECK_LT(left_node, num_left_nodes_); 00290 ArcIndex matching_arc = GetAssignmentArc(left_node); 00291 DCHECK_NE(GraphType::kNilArc, matching_arc); 00292 return Head(matching_arc); 00293 } 00294 00295 string StatsString() const { 00296 return total_stats_.StatsString(); 00297 } 00298 00299 class BipartiteLeftNodeIterator { 00300 public: 00301 BipartiteLeftNodeIterator(const GraphType& graph, NodeIndex num_left_nodes) 00302 : num_left_nodes_(num_left_nodes), 00303 node_iterator_(graph) { } 00304 00305 explicit BipartiteLeftNodeIterator(const LinearSumAssignment& assignment) 00306 : num_left_nodes_(assignment.NumLeftNodes()), 00307 node_iterator_(assignment.Graph()) { } 00308 00309 NodeIndex Index() const { return node_iterator_.Index(); } 00310 00311 bool Ok() const { 00312 return node_iterator_.Ok() && (node_iterator_.Index() < num_left_nodes_); 00313 } 00314 00315 void Next() { node_iterator_.Next(); } 00316 00317 private: 00318 const NodeIndex num_left_nodes_; 00319 typename GraphType::NodeIterator node_iterator_; 00320 }; 00321 00322 private: 00323 struct Stats { 00324 Stats() 00325 : pushes_(0), 00326 double_pushes_(0), 00327 relabelings_(0), 00328 refinements_(0) { } 00329 void Clear() { 00330 pushes_ = 0; 00331 double_pushes_ = 0; 00332 relabelings_ = 0; 00333 refinements_ = 0; 00334 } 00335 void Add(const Stats& that) { 00336 pushes_ += that.pushes_; 00337 double_pushes_ += that.double_pushes_; 00338 relabelings_ += that.relabelings_; 00339 refinements_ += that.refinements_; 00340 } 00341 string StatsString() const { 00342 return StringPrintf("%lld refinements; %lld relabelings; " 00343 "%lld double pushes; %lld pushes", 00344 refinements_, 00345 relabelings_, 00346 double_pushes_, 00347 pushes_); 00348 } 00349 int64 pushes_; 00350 int64 double_pushes_; 00351 int64 relabelings_; 00352 int64 refinements_; 00353 }; 00354 00355 #ifndef SWIG 00356 class ActiveNodeContainerInterface { 00357 public: 00358 virtual ~ActiveNodeContainerInterface() {} 00359 virtual bool Empty() const = 0; 00360 virtual void Add(NodeIndex node) = 0; 00361 virtual NodeIndex Get() = 0; 00362 }; 00363 00364 class ActiveNodeStack : public ActiveNodeContainerInterface { 00365 public: 00366 virtual ~ActiveNodeStack() {} 00367 00368 virtual bool Empty() const { 00369 return v_.empty(); 00370 } 00371 00372 virtual void Add(NodeIndex node) { 00373 v_.push_back(node); 00374 } 00375 00376 virtual NodeIndex Get() { 00377 DCHECK(!Empty()); 00378 NodeIndex result = v_.back(); 00379 v_.pop_back(); 00380 return result; 00381 } 00382 00383 private: 00384 std::vector<NodeIndex> v_; 00385 }; 00386 00387 class ActiveNodeQueue : public ActiveNodeContainerInterface { 00388 public: 00389 virtual ~ActiveNodeQueue() {} 00390 00391 virtual bool Empty() const { 00392 return q_.empty(); 00393 } 00394 00395 virtual void Add(NodeIndex node) { 00396 q_.push_front(node); 00397 } 00398 00399 virtual NodeIndex Get() { 00400 DCHECK(!Empty()); 00401 NodeIndex result= q_.back(); 00402 q_.pop_back(); 00403 return result; 00404 } 00405 00406 private: 00407 std::deque<NodeIndex> q_; 00408 }; 00409 #endif 00410 00411 // Type definition for a pair 00412 // (arc index, reduced cost gap) 00413 // giving the arc along which we will push from a given left-side 00414 // node and the gap between that arc's partial reduced cost and the 00415 // reduced cost of the next-best (necessarily residual) arc out of 00416 // the node. This information helps us efficiently relabel 00417 // right-side nodes during DoublePush operations. 00418 typedef std::pair<ArcIndex, CostValue> ImplicitPriceSummary; 00419 00420 // Returns true if and only if the current pseudoflow is 00421 // epsilon-optimal. To be used in a DCHECK. 00422 bool EpsilonOptimal() const; 00423 00424 // Checks that all nodes are matched. 00425 // To be used in a DCHECK. 00426 bool AllMatched() const; 00427 00428 // Calculates the implicit price of the given node. 00429 // Only for debugging, for use in EpsilonOptimal(). 00430 inline CostValue ImplicitPrice(NodeIndex left_node) const; 00431 00432 // For use by DoublePush() 00433 inline ImplicitPriceSummary BestArcAndGap(NodeIndex left_node) const; 00434 00435 // Accumulates stats between iterations and reports them if the 00436 // verbosity level is high enough. 00437 void ReportAndAccumulateStats(); 00438 00439 // Utility function to compute the next error parameter value. This 00440 // is used to ensure that the same sequence of error parameter 00441 // values is used for computation of price bounds as is used for 00442 // computing the optimum assignment. 00443 CostValue NewEpsilon(CostValue current_epsilon) const; 00444 00445 // Advances internal state to prepare for the next scaling 00446 // iteration. Returns false if infeasibility is detected, true 00447 // otherwise. 00448 bool UpdateEpsilon(); 00449 00450 // Indicates whether the given left_node has positive excess. Called 00451 // only for nodes on the left side. 00452 inline bool IsActive(NodeIndex left_node) const; 00453 00454 // Indicates whether the given node has nonzero excess. The idea 00455 // here is the same as the IsActive method above, but that method 00456 // contains a safety DCHECK() that its argument is a left-side node, 00457 // while this method is usable for any node. 00458 // To be used in a DCHECK. 00459 inline bool IsActiveForDebugging(NodeIndex node) const; 00460 00461 // Performs the push/relabel work for one scaling iteration. 00462 bool Refine(); 00463 00464 // Puts all left-side nodes in the active set in preparation for the 00465 // first scaling iteration. 00466 void InitializeActiveNodeContainer(); 00467 00468 // Saturates all negative-reduced-cost arcs at the beginning of each 00469 // scaling iteration. Note that according to the asymmetric 00470 // definition of admissibility, this action is different from 00471 // saturating all admissible arcs (which we never do). All negative 00472 // arcs are admissible, but not all admissible arcs are negative. It 00473 // is alwsys enough to saturate only the negative ones. 00474 void SaturateNegativeArcs(); 00475 00476 // Performs an optimized sequence of pushing a unit of excess out of 00477 // the left-side node v and back to another left-side node if no 00478 // deficit is cancelled with the first push. 00479 bool DoublePush(NodeIndex source); 00480 00481 // Returns the partial reduced cost of the given arc. 00482 inline CostValue PartialReducedCost(ArcIndex arc) const { 00483 return scaled_arc_cost_[arc] - price_[Head(arc)]; 00484 } 00485 00486 // The graph underlying the problem definition we are given. Not 00487 // const because we add arcs to the graph via our SetArcCost() 00488 // method. 00489 const GraphType& graph_; 00490 00491 // The number of nodes on the left side of the graph we are given. 00492 NodeIndex num_left_nodes_; 00493 00494 // A flag indicating that an optimal perfect matching has been computed. 00495 bool success_; 00496 00497 // The value by which we multiply all the arc costs we are given in 00498 // order to be able to use integer arithmetic in all our 00499 // computations. In order to establish optimality of the final 00500 // matching we compute, we need that 00501 // (cost_scaling_factor_ / kMinEpsilon) > graph_.num_nodes(). 00502 const CostValue cost_scaling_factor_; 00503 00504 // Scaling divisor. 00505 CostValue alpha_; 00506 00507 // Minimum value of epsilon. When a flow is epsilon-optimal for 00508 // epsilon == kMinEpsilon, the flow is optimal. 00509 static const CostValue kMinEpsilon; 00510 00511 // Current value of epsilon, the cost scaling parameter. 00512 CostValue epsilon_; 00513 00514 // The following two data members, price_lower_bound_ and 00515 // slack_relabeling_price_, have to do with bounds on the amount by 00516 // which node prices can change during execution of the algorithm. 00517 // We need some detailed discussion of this topic because we violate 00518 // several simplifying assumptions typically made in the theoretical 00519 // literature. In particular, we use integer arithmetic, we use a 00520 // reduction to the transportation problem rather than min-cost 00521 // circulation, we provide detection of infeasible problems rather 00522 // than assume feasibility, we detect when our computations might 00523 // exceed the range of representable cost values, and we use the 00524 // double-push heuristic which relabels nodes that do not have 00525 // excess. 00526 // 00527 // In the following discussion, we prove the following propositions: 00528 // Proposition 1. [Fidelity of arithmetic precision guarantee] If 00529 // FinalizeSetup() returns true, no arithmetic 00530 // overflow occurs during ComputeAssignment(). 00531 // Proposition 2. [Fidelity of feasibility detection] If no 00532 // arithmetic overflow occurs during 00533 // ComputeAssignment(), the return value of 00534 // ComputeAssignment() faithfully indicates whether 00535 // the given problem is feasible. 00536 // 00537 // We begin with some general discussion. 00538 // 00539 // The ideas used to prove our two propositions are essentially 00540 // those that appear in [Goldberg and Tarjan], but several details 00541 // are different: [Goldberg and Tarjan] assumes a feasible problem, 00542 // uses a symmetric notion of epsilon-optimality, considers only 00543 // nodes with excess eligible for relabeling, and does not treat the 00544 // question of arithmetic overflow. This implementation, on the 00545 // other hand, detects and reports infeasible problems, uses 00546 // asymmetric epsilon-optimality, relabels nodes with no excess in 00547 // the course of the double-push operation, and gives a reasonably 00548 // tight guarantee of arithmetic precision. No fundamentally new 00549 // ideas are involved, but the details are a bit tricky so they are 00550 // explained here. 00551 // 00552 // We have two intertwined needs that lead us to compute bounds on 00553 // the prices nodes can have during the assignment computation, on 00554 // the assumption that the given problem is feasible: 00555 // 1. Infeasibility detection: Infeasibility is detected by 00556 // observing that some node's price has been reduced too much by 00557 // relabeling operations (see [Goldberg and Tarjan] for the 00558 // argument -- duplicated in modified form below -- bounding the 00559 // running time of the push/relabel min-cost flow algorithm for 00560 // feasible problems); and 00561 // 2. Aggressively relabeling nodes and arcs whose matching is 00562 // forced: When a left-side node is incident to only one arc a, 00563 // any feasible solution must include a, and reducing the price 00564 // of Head(a) by any nonnegative amount preserves epsilon- 00565 // optimality. Because of this freedom, we'll call this sort of 00566 // relabeling (i.e., a relabeling of a right-side node that is 00567 // the only neighbor of the left-side node to which it has been 00568 // matched in the present double-push operation) a "slack" 00569 // relabeling. Relabelings that are not slack relabelings are 00570 // called "confined" relabelings. By relabeling Head(a) to have 00571 // p(Head(a))=-infinity, we could guarantee that a never becomes 00572 // unmatched during the current iteration, and this would prevent 00573 // our wasting time repeatedly unmatching and rematching a. But 00574 // there are some details we need to handle: 00575 // a. The CostValue type cannot represent -infinity; 00576 // b. Low node prices are precisely the signal we use to detect 00577 // infeasibility (see (1)), so we must be careful not to 00578 // falsely conclude that the problem is infeasible as a result 00579 // of the low price we gave Head(a); and 00580 // c. We need to indicate accurately to the client when our best 00581 // understanding indicates that we can't rule out arithmetic 00582 // overflow in our calculations. Most importantly, if we don't 00583 // warn the client, we must be certain to avoid overflow. This 00584 // means our slack relabelings must not be so aggressive as to 00585 // create the possibility of unforeseen overflow. Although we 00586 // will not achieve this in practice, slack relabelings would 00587 // ideally not introduce overflow unless overflow was 00588 // inevitable were even the smallest reasonable price change 00589 // (== epsilon) used for slack relabelings. 00590 // Using the analysis below, we choose a finite amount of price 00591 // change for slack relabelings aggressive enough that we don't 00592 // waste time doing repeated slack relabelings in a single 00593 // iteration, yet modest enough that we keep a good handle on 00594 // arithmetic precision and our ability to detect infeasible 00595 // problems. 00596 // 00597 // To provide faithful detection of infeasibility, a dependable 00598 // guarantee of arithmetic precision whenever possible, and good 00599 // performance by aggressively relabeling nodes whose matching is 00600 // forced, we exploit these facts: 00601 // 1. Beyond the first iteration, infeasibility detection isn't needed 00602 // because a problem is feasible in some iteration if and only if 00603 // it's feasible in all others. Therefore we are free to use an 00604 // infeasibility detection mechanism that might work in just one 00605 // iteration and switch it off in all other iterations. 00606 // 2. When we do a slack relabeling, we must choose the amount of 00607 // price reduction to use. We choose an amount large enough to 00608 // guarantee putting the node's matching to rest, yet (although 00609 // we don't bother to prove this explicitly) small enough that 00610 // the node's price obeys the overall lower bound that holds if 00611 // the slack relabeling amount is small. 00612 // 00613 // We will establish Propositions (1) and (2) above according to the 00614 // following steps: 00615 // First, we prove Lemma 1, which is a modified form of lemma 5.8 of 00616 // [Goldberg and Tarjan] giving a bound on the difference in price 00617 // between the end nodes of certain paths in the residual graph. 00618 // Second, we prove Lemma 2, which is technical lemma to establish 00619 // reachability of certain "anchor" nodes in the residual graph from 00620 // any node where a relabeling takes place. 00621 // Third, we apply the first two lemmas to prove Lemma 3 and Lemma 00622 // 4, which give two similar bounds that hold whenever the given 00623 // problem is feasible: (for feasibility detection) a bound on the 00624 // price of any node we relabel during any iteration (and the first 00625 // iteration in particular), and (for arithmetic precision) a bound 00626 // on the price of any node we relabel during the entire algorithm. 00627 // 00628 // Finally, we note that if the whole-algorithm price bound can be 00629 // represented precisely by the CostValue type, arithmetic overflow 00630 // cannot occur (establishing Proposition 1), and assuming no 00631 // overflow occurs during the first iteration, any violation of the 00632 // first-iteration price bound establishes infeasibility 00633 // (Proposition 2). 00634 // 00635 // The statement of Lemma 1 is perhaps easier to understand when the 00636 // reader knows how it will be used. To wit: In this lemma, f' and 00637 // e_0 are the flow and error parameter (epsilon) at the beginning 00638 // of the current iteration, while f and e_1 are the current 00639 // pseudoflow and error parameter when a relabeling of interest 00640 // occurs. Without loss of generality, c is the reduced cost 00641 // function at the beginning of the current iteration and p is the 00642 // change in prices that has taken place in the current iteration. 00643 // 00644 // Lemma 1 (a variant of lemma 5.8 from [Goldberg and Tarjan]): Let 00645 // f be a pseudoflow and let f' be a flow. Suppose P is a simple 00646 // path from right-side node v to right-side node w such that P is 00647 // residual with respect to f and reverse(P) is residual with 00648 // respect to f'. Further, suppose c is an arc cost function with 00649 // respect to which f' is e_0-optimal with the zero price function 00650 // and p is a price function with respect to which f is e_1-optimal 00651 // with respect to p. Then 00652 // p(v) - p(w) >= -(e_0 + e_1) * (n-2)/2. (***) 00653 // 00654 // Proof: We have c_p(P) = p(v) + c(P) - p(w) and hence 00655 // p(v) - p(w) = c_p(P) - c(P). 00656 // So we seek a bound on c_p(P) - c(P). 00657 // p(v) = c_p(P) - c(P). 00658 // Let arc a lie on P, which implies that a is residual with respect 00659 // to f and reverse(a) is residual with respect to f'. 00660 // Case 1: a is a forward arc. Then by e_1-optimality of f with 00661 // respect to p, c_p(a) >= 0 and reverse(a) is residual with 00662 // respect to f'. By e_0-optimality of f', c(a) <= e_0. So 00663 // c_p(a) - c(a) >= -e_0. 00664 // Case 2: a is a reverse arc. Then by e_1-optimality of f with 00665 // respect to p, c_p(a) >= -e_1 and reverse(a) is residual 00666 // with respect to f'. By e_0-optimality of f', c(a) <= 0. 00667 // So 00668 // c_p(a) - c(a) >= -e_1. 00669 // We assumed v and w are both right-side nodes, so there are at 00670 // most n - 2 arcs on the path P, of which at most (n-2)/2 are 00671 // forward arcs and at most (n-2)/2 are reverse arcs, so 00672 // p(v) - p(w) = c_p(P) - c(P) 00673 // >= -(e_0 + e_1) * (n-2)/2. (***) 00674 // 00675 // Some of the rest of our argument is given as a sketch, omitting 00676 // several details. Also elided here are some minor technical issues 00677 // related to the first iteration, inasmuch as our arguments assume 00678 // on the surface a "previous iteration" that doesn't exist in that 00679 // case. The issues are not substantial, just a bit messy. 00680 // 00681 // Lemma 2 is analogous to lemma 5.7 of [Goldberg and Tarjan], where 00682 // they have only relabelings that take place at nodes with excess 00683 // while we have only relabelings that take place as part of the 00684 // double-push operation at nodes without excess. 00685 // 00686 // Lemma 2: When a right-side node v is relabeled by our 00687 // implementation, either the problem is infeasible or there exists 00688 // a node w such that 00689 // A. w is reachable from v along some simple residual path P where 00690 // reverse(P) was residual at the beginning of the current 00691 // iteration; and 00692 // B. at least one of the following holds: 00693 // 1. when w was last relabeled, there existed a path P' from w 00694 // to a node with deficit in the residual graph where 00695 // reverse(P') was residual at the beginning of the current 00696 // iteration; or 00697 // 2. when w was last relabeled, it was a slack relabeling; 00698 // and 00699 // C. at least one of the following holds: 00700 // 1. w will not be relabeled again in this iteration; or 00701 // 2. v == w. 00702 // 00703 // The proof of Lemma 2 is somewhat messy and is omitted for 00704 // expedience. 00705 // 00706 // Lemma 1 bounds the price change during an iteration for any node 00707 // relabeled when a deficit is residually reachable from that node, 00708 // since a node w with deficit is not relabeled, hence p(w) = 0 in 00709 // the Lemma 1 bound. Let the bound from Lemma 1 with p(w) = 0 be 00710 // called B(e_0, e_1), and let us say that when a slack relabeling 00711 // of a node v occurs, we will set the price of v to B(e_0, e_1) 00712 // such that v tightly satisfies the bound of Lemma 1. Explicitly, 00713 // we define 00714 // B(e_0, e_1) = -(e_0 + e_1) * (n-2)/2. 00715 // 00716 // From Lemma 1 and Lemma 2, and taking into account our knowledge 00717 // of the slack relabeling amount, we have Lemma 3. 00718 // 00719 // Lemma 3: During any iteration, if the given problem is feasible 00720 // the price of any node is reduced by less than 00721 // 2 * B(e_0, e_1) = -(e_0 + e_1) * (n-2). 00722 // 00723 // Proof: Straightforward, omitted for expedience. 00724 // 00725 // In the case where e_0 = e_1 * alpha, we can express the bound 00726 // just in terms of e_1, the current iteration's value of epsilon_: 00727 // B(e_1) = B(e_1 * alpha, e_1) = -(1 + alpha) * e_1 * (n-2)/2, 00728 // so we have that p(v) is reduced by less than 2 * B(e_1). 00729 // 00730 // Because we use truncating division to compute each iteration's error 00731 // parameter from that of the previous iteration, it isn't exactly 00732 // the case that e_0 = e_1 * alpha as we just assumed. To patch this 00733 // up, we can use the observation that 00734 // e_1 = floor(e_0 / alpha), 00735 // which implies 00736 // -e_0 > -(e_1 + 1) * alpha 00737 // to rewrite from (***): 00738 // p(v) > 2 * B(e_0, e_1) > 2 * B((e_1 + 1) * alpha, e_1) 00739 // = 2 * -((e_1 + 1) * alpha + e_1) * (n-2)/2 00740 // = 2 * -(1 + alpha) * e_1 * (n-2)/2 - alpha * (n-2) 00741 // = 2 * B(e_1) - alpha * (n-2) 00742 // = -((1 + alpha) * e_1 + alpha) * (n-2). 00743 // 00744 // We sum up the bounds for all the iterations to get Lemma 4: 00745 // 00746 // Lemma 4: If the given problem is feasible, after k iterations the 00747 // price of any node is always greater than 00748 // -((1 + alpha) * C + (k * alpha)) * (n-2) 00749 // 00750 // Proof: Suppose the price decrease of every node in the iteration 00751 // with epsilon_ == x is bounded by B(x) which is proportional to x 00752 // (not surpisingly, this will be the same function B() as 00753 // above). Assume for simplicity that C, the largest cost magnitude, 00754 // is a power of alpha. Then the price of each node, tallied across 00755 // all iterations is bounded 00756 // p(v) > 2 * B(C/alpha) + 2 * B(C/alpha^2) + ... + 2 * B(kMinEpsilon) 00757 // == 2 * B(C/alpha) * alpha / (alpha - 1) 00758 // == 2 * B(C) / (alpha - 1). 00759 // As above, this needs some patching up to handle the fact that we 00760 // use truncating arithmetic. We saw that each iteration effectively 00761 // reduces the price bound by alpha * (n-2), hence if there are k 00762 // iterations, the bound is 00763 // p(v) > 2 * B(C) / (alpha - 1) - k * alpha * (n-2) 00764 // = -(1 + alpha) * C * (n-2) / (alpha - 1) - k * alpha * (n-2) 00765 // = (n-2) * (C * (1 + alpha) / (1 - alpha) - k * alpha). 00766 // 00767 // The bound of lemma 4 can be used to warn for possible overflow of 00768 // arithmetic precision. But because it involves the number of 00769 // iterations, k, we might as well count through the iterations 00770 // simply adding up the bounds given by Lemma 3 to get a tighter 00771 // result. This is what the implementation does. 00772 00773 // A lower bound on the price of any node at any time throughout the 00774 // computation. A price below this level proves infeasibility; this 00775 // value is used for feasibility detection. We use this value also 00776 // to rule out the possibility of arithmetic overflow or warn the 00777 // client that we have not been able to rule out that possibility. 00778 // 00779 // We can use the value implied by Lemma 4 here, but note that that 00780 // value includes k, the number of iterations. It's plenty fast if 00781 // we count through the iterations to compute that value, but if 00782 // we're going to count through the iterations, we might as well use 00783 // the two-parameter bound from Lemma 3, summing up as we go. This 00784 // gives us a tighter bound and more comprehensible code. 00785 // 00786 // While computing this bound, if we find the value justified by the 00787 // theory lies outside the representable range of CostValue, we 00788 // conclude that the given arc costs have magnitudes so large that 00789 // we cannot guarantee our calculations don't overflow. If the value 00790 // justified by the theory lies inside the representable range of 00791 // CostValue, we commit that our calculation will not overflow. This 00792 // commitment means we need to be careful with the amount by which 00793 // we relabel right-side nodes that are incident to any node with 00794 // only one neighbor. 00795 CostValue price_lower_bound_; 00796 00797 // A bound on the amount by which a node's price can be reduced 00798 // during the current iteration, used only for slack 00799 // relabelings. Where epsilon is the first iteration's error 00800 // parameter and C is the largest magnitude of an arc cost, we set 00801 // slack_relabeling_price_ = -B(C, epsilon) 00802 // = (C + epsilon) * (n-2)/2. 00803 // 00804 // We could use slack_relabeling_price_ for feasibility detection 00805 // but the feasibility threshold is double the slack relabeling 00806 // amount and we judge it not to be worth having to multiply by two 00807 // gratuitously to check feasibility in each double push 00808 // operation. Instead we settle for feasibility detection using 00809 // price_lower_bound_ instead, which is somewhat slower in the 00810 // infeasible case because more relabelings will be required for 00811 // some node price to attain the looser bound. 00812 CostValue slack_relabeling_price_; 00813 00814 // Computes the value of the bound on price reduction for an 00815 // iteration, given the old and new values of epsilon_. Because the 00816 // expression computed here is used in at least one place where we 00817 // want an additional factor in the denominator, we take that factor 00818 // as an argument. If extra_divisor == 1, this function computes of 00819 // the function B() discussed above. 00820 // 00821 // Avoids overflow in computing the bound, and sets *in_range = 00822 // false if the value of the bound doesn't fit in CostValue. 00823 inline CostValue PriceChangeBound(CostValue old_epsilon, 00824 CostValue new_epsilon, 00825 bool* in_range) const { 00826 const CostValue n = graph_.num_nodes(); 00827 // We work in double-precision floating point to determine whether 00828 // we'll overflow the integral CostValue type's range of 00829 // representation. Switching between integer and double is a 00830 // rather expensive operation, but we do this only twice per 00831 // scaling iteration, so we can afford it rather than resort to 00832 // complex and subtle tricks within the bounds of integer 00833 // arithmetic. 00834 // 00835 // You will want to read the comments above about 00836 // price_lower_bound_ and slack_relabeling_price_, and have a 00837 // pencil handy. :-) 00838 const double result = 00839 static_cast<double>(std::max<CostValue>(0, n / 2 - 1)) * 00840 static_cast<double>(old_epsilon + new_epsilon); 00841 const double limit = 00842 static_cast<double>(std::numeric_limits<CostValue>::max()); 00843 if (result > limit) { 00844 // Our integer computations could overflow. 00845 if (in_range != NULL) *in_range = false; 00846 return std::numeric_limits<CostValue>::max(); 00847 } else { 00848 // Don't touch *in_range; other computations could already have 00849 // set it to false and we don't want to overwrite that result. 00850 return static_cast<CostValue>(result); 00851 } 00852 } 00853 00854 // A scaled record of the largest arc-cost magnitude we've been 00855 // given during problem setup. This is used to set the initial value 00856 // of epsilon_, which in turn is used not only as the error 00857 // parameter but also to determine whether we risk arithmetic 00858 // overflow during the algorithm. 00859 // 00860 // Note: Our treatment of arithmetic overflow assumes the following 00861 // property of CostValue: 00862 // -std::numeric_limits<CostValue>::max() is a representable 00863 // CostValue. 00864 // That property is satisfied if CostValue uses a two's-complement 00865 // representation. 00866 CostValue largest_scaled_cost_magnitude_; 00867 00868 // The total excess in the graph. Given our asymmetric definition of 00869 // epsilon-optimality and our use of the double-push operation, this 00870 // equals the number of unmatched left-side nodes. 00871 NodeIndex total_excess_; 00872 00873 // Indexed by node index, the price_ values are maintained only for 00874 // right-side nodes. 00875 CostArray price_; 00876 00877 // Indexed by left-side node index, the matched_arc_ array gives the 00878 // arc index of the arc matching any given left-side node, or 00879 // GraphType::kNilArc if the node is unmatched. 00880 ArcIndexArray matched_arc_; 00881 00882 // Indexed by right-side node index, the matched_node_ array gives 00883 // the node index of the left-side node matching any given 00884 // right-side node, or GraphType::kNilNode if the right-side node is 00885 // unmatched. 00886 NodeIndexArray matched_node_; 00887 00888 // The array of arc costs as given in the problem definition, except 00889 // that they are scaled up by the number of nodes in the graph so we 00890 // can use integer arithmetic throughout. 00891 CostArray scaled_arc_cost_; 00892 00893 // The container of active nodes (i.e., unmatched nodes). This can 00894 // be switched easily between ActiveNodeStack and ActiveNodeQueue 00895 // for experimentation. 00896 scoped_ptr<ActiveNodeContainerInterface> active_nodes_; 00897 00898 // Statistics giving the overall numbers of various operations the 00899 // algorithm performs. 00900 Stats total_stats_; 00901 00902 // Statistics giving the numbers of various operations the algorithm 00903 // has performed in the current iteration. 00904 Stats iteration_stats_; 00905 00906 DISALLOW_COPY_AND_ASSIGN(LinearSumAssignment); 00907 }; 00908 00909 // Implementation of out-of-line LinearSumAssignment template member 00910 // functions. 00911 00912 template <typename GraphType> 00913 const CostValue LinearSumAssignment<GraphType>::kMinEpsilon = 1; 00914 00915 template <typename GraphType> 00916 LinearSumAssignment<GraphType>::LinearSumAssignment( 00917 const GraphType& graph, NodeIndex num_left_nodes) 00918 : graph_(graph), 00919 num_left_nodes_(num_left_nodes), 00920 success_(false), 00921 cost_scaling_factor_(1 + (graph.max_num_nodes() / 2)), 00922 alpha_(FLAGS_assignment_alpha), 00923 epsilon_(0), 00924 price_lower_bound_(0), 00925 slack_relabeling_price_(0), 00926 largest_scaled_cost_magnitude_(0), 00927 total_excess_(0), 00928 price_(num_left_nodes + GraphType::kFirstNode, 00929 graph.max_end_node_index() - 1), 00930 matched_arc_(GraphType::kFirstNode, num_left_nodes - 1), 00931 matched_node_(num_left_nodes, graph.max_end_node_index() - 1), 00932 scaled_arc_cost_(GraphType::kFirstArc, graph.max_end_arc_index() - 1), 00933 active_nodes_( 00934 FLAGS_assignment_stack_order ? 00935 static_cast<ActiveNodeContainerInterface*>(new ActiveNodeStack()) : 00936 static_cast<ActiveNodeContainerInterface*>(new ActiveNodeQueue())) { } 00937 00938 template <typename GraphType> 00939 void LinearSumAssignment<GraphType>::SetArcCost(ArcIndex arc, CostValue cost) { 00940 DCHECK(graph_.CheckArcValidity(arc)); 00941 NodeIndex head = Head(arc); 00942 DCHECK_LE(num_left_nodes_, head); 00943 cost *= cost_scaling_factor_; 00944 const CostValue cost_magnitude = std::abs(cost); 00945 largest_scaled_cost_magnitude_ = std::max(largest_scaled_cost_magnitude_, 00946 cost_magnitude); 00947 scaled_arc_cost_.Set(arc, cost); 00948 } 00949 00950 template <typename ArcIndexType> 00951 class CostValueCycleHandler 00952 : public PermutationCycleHandler<ArcIndexType> { 00953 public: 00954 explicit CostValueCycleHandler(CostArray* cost) 00955 : temp_(0), 00956 cost_(cost) { } 00957 00958 virtual void SetTempFromIndex(ArcIndexType source) { 00959 temp_ = cost_->Value(source); 00960 } 00961 00962 virtual void SetIndexFromIndex(ArcIndexType source, 00963 ArcIndexType destination) const { 00964 cost_->Set(destination, cost_->Value(source)); 00965 } 00966 00967 virtual void SetIndexFromTemp(ArcIndexType destination) const { 00968 cost_->Set(destination, temp_); 00969 } 00970 00971 virtual ~CostValueCycleHandler() { } 00972 00973 private: 00974 CostValue temp_; 00975 00976 CostArray* cost_; 00977 00978 DISALLOW_COPY_AND_ASSIGN(CostValueCycleHandler); 00979 }; 00980 00981 // Logically this class should be defined inside OptimizeGraphLayout, 00982 // but compilation fails if we do that because C++98 doesn't allow 00983 // instantiation of member templates with function-scoped types as 00984 // template parameters, which in turn is because those function-scoped 00985 // types lack linkage. 00986 template <typename GraphType> class ArcIndexOrderingByTailNode { 00987 public: 00988 explicit ArcIndexOrderingByTailNode(const GraphType& graph) 00989 : graph_(graph) { } 00990 00991 // Says ArcIndex a is less than ArcIndex b if arc a's tail is less 00992 // than arc b's tail. If their tails are equal, orders according to 00993 // heads. 00994 bool operator()(ArcIndex a, ArcIndex b) const { 00995 return ((graph_.Tail(a) < graph_.Tail(b)) || 00996 ((graph_.Tail(a) == graph_.Tail(b)) && 00997 (graph_.Head(a) < graph_.Head(b)))); 00998 } 00999 01000 private: 01001 const GraphType& graph_; 01002 01003 // Copy and assign are allowed; they have to be for STL to work 01004 // with this functor, although it seems like a bug for STL to be 01005 // written that way. 01006 }; 01007 01008 template <typename GraphType> 01009 void LinearSumAssignment<GraphType>::OptimizeGraphLayout(GraphType* graph) { 01010 // The graph argument is only to give us a non-const-qualified 01011 // handle on the graph we already have. Any different graph is 01012 // nonsense. 01013 DCHECK_EQ(&graph_, graph); 01014 const ArcIndexOrderingByTailNode<GraphType> compare(graph_); 01015 CostValueCycleHandler<typename GraphType::ArcIndex> 01016 cycle_handler(&scaled_arc_cost_); 01017 TailArrayManager<GraphType> tail_array_manager(graph); 01018 tail_array_manager.BuildTailArrayFromAdjacencyListsIfForwardGraph(); 01019 graph->GroupForwardArcsByFunctor(compare, &cycle_handler); 01020 tail_array_manager.ReleaseTailArrayIfForwardGraph(); 01021 } 01022 01023 template <typename GraphType> 01024 CostValue LinearSumAssignment<GraphType>::NewEpsilon( 01025 const CostValue current_epsilon) const { 01026 return std::max(current_epsilon / alpha_, kMinEpsilon); 01027 } 01028 01029 template <typename GraphType> 01030 bool LinearSumAssignment<GraphType>::UpdateEpsilon() { 01031 CostValue new_epsilon = NewEpsilon(epsilon_); 01032 slack_relabeling_price_ = PriceChangeBound(epsilon_, new_epsilon, NULL); 01033 epsilon_ = new_epsilon; 01034 VLOG(3) << "Updated: epsilon_ == " << epsilon_; 01035 VLOG(4) << "slack_relabeling_price_ == " << slack_relabeling_price_; 01036 DCHECK_GT(slack_relabeling_price_, 0); 01037 // For today we always return true; in the future updating epsilon 01038 // in sophisticated ways could conceivably detect infeasibility 01039 // before the first iteration of Refine(). 01040 return true; 01041 } 01042 01043 // For production code that checks whether a left-side node is active. 01044 template <typename GraphType> 01045 inline bool LinearSumAssignment<GraphType>::IsActive( 01046 NodeIndex left_node) const { 01047 DCHECK_LT(left_node, num_left_nodes_); 01048 return matched_arc_[left_node] == GraphType::kNilArc; 01049 } 01050 01051 // Only for debugging. Separate from the production IsActive() method 01052 // so that method can assert that its argument is a left-side node, 01053 // while for debugging we need to be able to test any node. 01054 template <typename GraphType> 01055 inline bool LinearSumAssignment<GraphType>::IsActiveForDebugging( 01056 NodeIndex node) const { 01057 if (node < num_left_nodes_) { 01058 return IsActive(node); 01059 } else { 01060 return matched_node_[node] == GraphType::kNilNode; 01061 } 01062 } 01063 01064 template <typename GraphType> 01065 void LinearSumAssignment<GraphType>::InitializeActiveNodeContainer() { 01066 DCHECK(active_nodes_->Empty()); 01067 for (BipartiteLeftNodeIterator node_it(graph_, num_left_nodes_); 01068 node_it.Ok(); 01069 node_it.Next()) { 01070 const NodeIndex node = node_it.Index(); 01071 if (IsActive(node)) { 01072 active_nodes_->Add(node); 01073 } 01074 } 01075 } 01076 01077 // There exists a price function such that the admissible arcs at the 01078 // beginning of an iteration are exactly the reverse arcs of all 01079 // matching arcs. Saturating all admissible arcs with respect to that 01080 // price function therefore means simply unmatching every matched 01081 // node. 01082 // 01083 // In the future we will price out arcs, which will reduce the set of 01084 // nodes we unmatch here. If a matching arc is priced out, we will not 01085 // unmatch its endpoints since that element of the matching is 01086 // guaranteed not to change. 01087 template <typename GraphType> 01088 void LinearSumAssignment<GraphType>::SaturateNegativeArcs() { 01089 total_excess_ = 0; 01090 for (BipartiteLeftNodeIterator node_it(graph_, num_left_nodes_); 01091 node_it.Ok(); 01092 node_it.Next()) { 01093 const NodeIndex node = node_it.Index(); 01094 if (IsActive(node)) { 01095 // This can happen in the first iteration when nothing is 01096 // matched yet. 01097 total_excess_ += 1; 01098 } else { 01099 // We're about to create a unit of excess by unmatching these nodes. 01100 total_excess_ += 1; 01101 const NodeIndex mate = GetMate(node); 01102 matched_arc_.Set(node, GraphType::kNilArc); 01103 matched_node_.Set(mate, GraphType::kNilNode); 01104 } 01105 } 01106 } 01107 01108 // Returns true for success, false for infeasible. 01109 template <typename GraphType> 01110 bool LinearSumAssignment<GraphType>::DoublePush(NodeIndex source) { 01111 DCHECK_GT(num_left_nodes_, source); 01112 DCHECK(IsActive(source)); 01113 ImplicitPriceSummary summary = BestArcAndGap(source); 01114 const ArcIndex best_arc = summary.first; 01115 const CostValue gap = summary.second; 01116 // Now we have the best arc incident to source, i.e., the one with 01117 // minimum reduced cost. Match that arc, unmatching its head if 01118 // necessary. 01119 if (best_arc == GraphType::kNilArc) { 01120 return false; 01121 } 01122 const NodeIndex new_mate = Head(best_arc); 01123 const NodeIndex to_unmatch = matched_node_[new_mate]; 01124 if (to_unmatch != GraphType::kNilNode) { 01125 // Unmatch new_mate from its current mate, pushing the unit of 01126 // flow back to a node on the left side as a unit of excess. 01127 matched_arc_.Set(to_unmatch, GraphType::kNilArc); 01128 active_nodes_->Add(to_unmatch); 01129 // This counts as a double push. 01130 iteration_stats_.double_pushes_ += 1; 01131 } else { 01132 // We are about to increase the cardinality of the matching. 01133 total_excess_ -= 1; 01134 // This counts as a single push. 01135 iteration_stats_.pushes_ += 1; 01136 } 01137 matched_arc_.Set(source, best_arc); 01138 matched_node_.Set(new_mate, source); 01139 // Finally, relabel new_mate. 01140 iteration_stats_.relabelings_ += 1; 01141 CostValue new_price = price_[new_mate] - gap - epsilon_; 01142 price_.Set(new_mate, new_price); 01143 return new_price >= price_lower_bound_; 01144 } 01145 01146 template <typename GraphType> 01147 bool LinearSumAssignment<GraphType>::Refine() { 01148 SaturateNegativeArcs(); 01149 InitializeActiveNodeContainer(); 01150 while (total_excess_ > 0) { 01151 // Get an active node (i.e., one with excess == 1) and discharge 01152 // it using DoublePush. 01153 const NodeIndex node = active_nodes_->Get(); 01154 if (!DoublePush(node)) { 01155 // Infeasibility detected. 01156 return false; 01157 } 01158 } 01159 DCHECK(active_nodes_->Empty()); 01160 iteration_stats_.refinements_ += 1; 01161 return true; 01162 } 01163 01164 // Computes best_arc, the minimum reduced-cost arc incident to 01165 // left_node and admissibility_gap, the amount by which the reduced 01166 // cost of best_arc must be increased to make it equal in reduced cost 01167 // to another residual arc incident to left_node. 01168 // 01169 // Precondition: left_node is unmatched. This allows us to simplify 01170 // the code. The debug-only counterpart to this routine is 01171 // LinearSumAssignment::ImplicitPrice() and it does not assume this 01172 // precondition. 01173 // 01174 // This function is large enough that our suggestion that the compiler 01175 // inline it might be pointless. 01176 template <typename GraphType> 01177 inline typename LinearSumAssignment<GraphType>::ImplicitPriceSummary 01178 LinearSumAssignment<GraphType>::BestArcAndGap(NodeIndex left_node) const { 01179 DCHECK(IsActive(left_node)); 01180 DCHECK_GT(epsilon_, 0); 01181 typename GraphType::OutgoingArcIterator arc_it(graph_, left_node); 01182 ArcIndex best_arc = arc_it.Index(); 01183 CostValue min_partial_reduced_cost = PartialReducedCost(best_arc); 01184 // We choose second_min_partial_reduced_cost so that in the case of 01185 // the largest possible gap (which results from a left-side node 01186 // with only a single incident residual arc), the corresponding 01187 // right-side node will be relabeled by an amount that exactly 01188 // matches slack_relabeling_price_. 01189 CostValue second_min_partial_reduced_cost = 01190 min_partial_reduced_cost + slack_relabeling_price_ - epsilon_; 01191 for (arc_it.Next(); arc_it.Ok(); arc_it.Next()) { 01192 const ArcIndex arc = arc_it.Index(); 01193 const CostValue partial_reduced_cost = PartialReducedCost(arc); 01194 if (partial_reduced_cost < second_min_partial_reduced_cost) { 01195 if (partial_reduced_cost < min_partial_reduced_cost) { 01196 best_arc = arc; 01197 second_min_partial_reduced_cost = min_partial_reduced_cost; 01198 min_partial_reduced_cost = partial_reduced_cost; 01199 } else { 01200 second_min_partial_reduced_cost = partial_reduced_cost; 01201 } 01202 } 01203 } 01204 const CostValue gap = 01205 second_min_partial_reduced_cost - min_partial_reduced_cost; 01206 DCHECK_GE(gap, 0); 01207 return std::make_pair(best_arc, gap); 01208 } 01209 01210 // Only for debugging. 01211 template <typename GraphType> inline CostValue 01212 LinearSumAssignment<GraphType>::ImplicitPrice(NodeIndex left_node) const { 01213 DCHECK_GT(num_left_nodes_, left_node); 01214 DCHECK_GT(epsilon_, 0); 01215 typename GraphType::OutgoingArcIterator arc_it(graph_, left_node); 01216 // If the input problem is feasible, it is always the case that 01217 // arc_it.Ok(), i.e., that there is at least one arc incident to 01218 // left_node. 01219 DCHECK(arc_it.Ok()); 01220 ArcIndex best_arc = arc_it.Index(); 01221 if (best_arc == matched_arc_[left_node]) { 01222 arc_it.Next(); 01223 if (arc_it.Ok()) { 01224 best_arc = arc_it.Index(); 01225 } 01226 } 01227 CostValue min_partial_reduced_cost = PartialReducedCost(best_arc); 01228 if (!arc_it.Ok()) { 01229 // Only one arc is incident to left_node, and the node is 01230 // currently matched along that arc, which must be the case in any 01231 // feasible solution. Therefore we implicitly price this node so 01232 // low that we will never consider unmatching it. 01233 return -(min_partial_reduced_cost + slack_relabeling_price_); 01234 } 01235 for (arc_it.Next(); arc_it.Ok(); arc_it.Next()) { 01236 const ArcIndex arc = arc_it.Index(); 01237 if (arc != matched_arc_[left_node]) { 01238 const CostValue partial_reduced_cost = PartialReducedCost(arc); 01239 if (partial_reduced_cost < min_partial_reduced_cost) { 01240 min_partial_reduced_cost = partial_reduced_cost; 01241 } 01242 } 01243 } 01244 return -min_partial_reduced_cost; 01245 } 01246 01247 // Only for debugging. 01248 template <typename GraphType> 01249 bool LinearSumAssignment<GraphType>::AllMatched() const { 01250 for (typename GraphType::NodeIterator node_it(graph_); 01251 node_it.Ok(); 01252 node_it.Next()) { 01253 if (IsActiveForDebugging(node_it.Index())) { 01254 return false; 01255 } 01256 } 01257 return true; 01258 } 01259 01260 // Only for debugging. 01261 template <typename GraphType> 01262 bool LinearSumAssignment<GraphType>::EpsilonOptimal() const { 01263 for (BipartiteLeftNodeIterator node_it(graph_, num_left_nodes_); 01264 node_it.Ok(); 01265 node_it.Next()) { 01266 const NodeIndex left_node = node_it.Index(); 01267 // Get the implicit price of left_node and make sure the reduced 01268 // costs of left_node's incident arcs are in bounds. 01269 CostValue left_node_price = ImplicitPrice(left_node); 01270 for (typename GraphType::OutgoingArcIterator arc_it(graph_, left_node); 01271 arc_it.Ok(); 01272 arc_it.Next()) { 01273 const ArcIndex arc = arc_it.Index(); 01274 const CostValue reduced_cost = 01275 left_node_price + PartialReducedCost(arc); 01276 // Note the asymmetric definition of epsilon-optimality that we 01277 // use because it means we can saturate all admissible arcs in 01278 // the beginning of Refine() just by unmatching all matched 01279 // nodes. 01280 if (matched_arc_[left_node] == arc) { 01281 // The reverse arc is residual. Epsilon-optimality requires 01282 // that the reduced cost of the forward arc be at most 01283 // epsilon_. 01284 if (reduced_cost > epsilon_) { 01285 return false; 01286 } 01287 } else { 01288 // The forward arc is residual. Epsilon-optimality requires 01289 // that the reduced cost of the forward arc be at least zero. 01290 if (reduced_cost < 0) { 01291 return false; 01292 } 01293 } 01294 } 01295 } 01296 return true; 01297 } 01298 01299 template <typename GraphType> 01300 bool LinearSumAssignment<GraphType>::FinalizeSetup() { 01301 epsilon_ = largest_scaled_cost_magnitude_; 01302 VLOG(2) << "Largest given cost magnitude: " << 01303 largest_scaled_cost_magnitude_ / cost_scaling_factor_; 01304 // Initialize left-side node-indexed arrays. 01305 typename GraphType::NodeIterator node_it(graph_); 01306 for (; node_it.Ok(); node_it.Next()) { 01307 const NodeIndex node = node_it.Index(); 01308 if (node >= num_left_nodes_) { 01309 break; 01310 } 01311 matched_arc_.Set(node, GraphType::kNilArc); 01312 } 01313 // Initialize right-side node-indexed arrays. Example: prices are 01314 // stored only for right-side nodes. 01315 for (; node_it.Ok(); node_it.Next()) { 01316 const NodeIndex node = node_it.Index(); 01317 price_.Set(node, 0); 01318 matched_node_.Set(node, GraphType::kNilNode); 01319 } 01320 bool in_range = true; 01321 double double_price_lower_bound = 0.0; 01322 CostValue new_error_parameter; 01323 CostValue old_error_parameter = epsilon_; 01324 do { 01325 new_error_parameter = NewEpsilon(old_error_parameter); 01326 double_price_lower_bound -= 2.0 * PriceChangeBound(old_error_parameter, 01327 new_error_parameter, 01328 &in_range); 01329 old_error_parameter = new_error_parameter; 01330 } while (new_error_parameter != kMinEpsilon); 01331 const double limit = 01332 -static_cast<double>(std::numeric_limits<CostValue>::max()); 01333 if (double_price_lower_bound < limit) { 01334 in_range = false; 01335 price_lower_bound_ = -std::numeric_limits<CostValue>::max(); 01336 } else { 01337 price_lower_bound_ = static_cast<CostValue>(double_price_lower_bound); 01338 } 01339 VLOG(4) << "price_lower_bound_ == " << price_lower_bound_; 01340 DCHECK_LE(price_lower_bound_, 0); 01341 if (!in_range) { 01342 LOG(WARNING) << "Price change bound exceeds range of representable " 01343 << "costs; arithmetic overflow is not ruled out and " 01344 << "infeasibility might go undetected."; 01345 } 01346 return in_range; 01347 } 01348 01349 template <typename GraphType> 01350 void LinearSumAssignment<GraphType>::ReportAndAccumulateStats() { 01351 total_stats_.Add(iteration_stats_); 01352 VLOG(3) << "Iteration stats: " << iteration_stats_.StatsString(); 01353 iteration_stats_.Clear(); 01354 } 01355 01356 template <typename GraphType> 01357 bool LinearSumAssignment<GraphType>::ComputeAssignment() { 01358 // Note: FinalizeSetup() might have been called already by white-box 01359 // test code or by a client that wants to react to the possibility 01360 // of overflow before solving the given problem, but FinalizeSetup() 01361 // is idempotent and reasonably fast, so we call it unconditionally 01362 // here. 01363 FinalizeSetup(); 01364 bool ok = graph_.num_nodes() == 2 * num_left_nodes_; 01365 DCHECK(!ok || EpsilonOptimal()); 01366 while (ok && epsilon_ > kMinEpsilon) { 01367 ok &= UpdateEpsilon(); 01368 ok &= Refine(); 01369 ReportAndAccumulateStats(); 01370 DCHECK(!ok || EpsilonOptimal()); 01371 DCHECK(!ok || AllMatched()); 01372 } 01373 success_ = ok; 01374 VLOG(1) << "Overall stats: " << total_stats_.StatsString(); 01375 return ok; 01376 } 01377 01378 template <typename GraphType> 01379 CostValue LinearSumAssignment<GraphType>::GetCost() const { 01380 // It is illegal to call this method unless we successfully computed 01381 // an optimum assignment. 01382 DCHECK(success_); 01383 CostValue cost = 0; 01384 for (BipartiteLeftNodeIterator node_it(*this); 01385 node_it.Ok(); 01386 node_it.Next()) { 01387 cost += GetAssignmentCost(node_it.Index()); 01388 } 01389 return cost; 01390 } 01391 01392 } // namespace operations_research 01393 01394 #endif // OR_TOOLS_GRAPH_LINEAR_ASSIGNMENT_H_