1 /* -*- mode: C++; indent-tabs-mode: nil; -*- 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 2003-2013 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 * 9 * Permission to use, modify and distribute this software is granted 10 * provided that this copyright notice appears in all copies. For 11 * precise terms see the accompanying LICENSE file. 12 * 13 * This software is provided "AS IS" with no warranty of any kind, 14 * express or implied, and with no claim as to its suitability for any 15 * purpose. 16 * 17 */ 18 19 #ifndef LEMON_CYCLE_CANCELING_H 20 #define LEMON_CYCLE_CANCELING_H 21 22 /// \ingroup min_cost_flow_algs 23 /// \file 24 /// \brief Cycle-canceling algorithms for finding a minimum cost flow. 25 26 #include <vector> 27 #include <limits> 28 29 #include <lemon/core.h> 30 #include <lemon/maps.h> 31 #include <lemon/path.h> 32 #include <lemon/math.h> 33 #include <lemon/static_graph.h> 34 #include <lemon/adaptors.h> 35 #include <lemon/circulation.h> 36 #include <lemon/bellman_ford.h> 37 #include <lemon/howard_mmc.h> 38 #include <lemon/hartmann_orlin_mmc.h> 39 40 namespace lemon { 41 42 /// \addtogroup min_cost_flow_algs 43 /// @{ 44 45 /// \brief Implementation of cycle-canceling algorithms for 46 /// finding a \ref min_cost_flow "minimum cost flow". 47 /// 48 /// \ref CycleCanceling implements three different cycle-canceling 49 /// algorithms for finding a \ref min_cost_flow "minimum cost flow" 50 /// \cite amo93networkflows, \cite klein67primal, 51 /// \cite goldberg89cyclecanceling. 52 /// The most efficent one is the \ref CANCEL_AND_TIGHTEN 53 /// "Cancel-and-Tighten" algorithm, thus it is the default method. 54 /// It runs in strongly polynomial time \f$O(n^2 m^2 \log n)\f$, 55 /// but in practice, it is typically orders of magnitude slower than 56 /// the scaling algorithms and \ref NetworkSimplex. 57 /// (For more information, see \ref min_cost_flow_algs "the module page".) 58 /// 59 /// Most of the parameters of the problem (except for the digraph) 60 /// can be given using separate functions, and the algorithm can be 61 /// executed using the \ref run() function. If some parameters are not 62 /// specified, then default values will be used. 63 /// 64 /// \tparam GR The digraph type the algorithm runs on. 65 /// \tparam V The number type used for flow amounts, capacity bounds 66 /// and supply values in the algorithm. By default, it is \c int. 67 /// \tparam C The number type used for costs and potentials in the 68 /// algorithm. By default, it is the same as \c V. 69 /// 70 /// \warning Both \c V and \c C must be signed number types. 71 /// \warning All input data (capacities, supply values, and costs) must 72 /// be integer. 73 /// \warning This algorithm does not support negative costs for 74 /// arcs having infinite upper bound. 75 /// 76 /// \note For more information about the three available methods, 77 /// see \ref Method. 78 #ifdef DOXYGEN 79 template <typename GR, typename V, typename C> 80 #else 81 template <typename GR, typename V = int, typename C = V> 82 #endif 83 class CycleCanceling 84 { 85 public: 86 87 /// The type of the digraph 88 typedef GR Digraph; 89 /// The type of the flow amounts, capacity bounds and supply values 90 typedef V Value; 91 /// The type of the arc costs 92 typedef C Cost; 93 94 public: 95 96 /// \brief Problem type constants for the \c run() function. 97 /// 98 /// Enum type containing the problem type constants that can be 99 /// returned by the \ref run() function of the algorithm. 100 enum ProblemType { 101 /// The problem has no feasible solution (flow). 102 INFEASIBLE, 103 /// The problem has optimal solution (i.e. it is feasible and 104 /// bounded), and the algorithm has found optimal flow and node 105 /// potentials (primal and dual solutions). 106 OPTIMAL, 107 /// The digraph contains an arc of negative cost and infinite 108 /// upper bound. It means that the objective function is unbounded 109 /// on that arc, however, note that it could actually be bounded 110 /// over the feasible flows, but this algroithm cannot handle 111 /// these cases. 112 UNBOUNDED 113 }; 114 115 /// \brief Constants for selecting the used method. 116 /// 117 /// Enum type containing constants for selecting the used method 118 /// for the \ref run() function. 119 /// 120 /// \ref CycleCanceling provides three different cycle-canceling 121 /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel-and-Tighten" 122 /// is used, which is by far the most efficient and the most robust. 123 /// However, the other methods can be selected using the \ref run() 124 /// function with the proper parameter. 125 enum Method { 126 /// A simple cycle-canceling method, which uses the 127 /// \ref BellmanFord "Bellman-Ford" algorithm for detecting negative 128 /// cycles in the residual network. 129 /// The number of Bellman-Ford iterations is bounded by a successively 130 /// increased limit. 131 SIMPLE_CYCLE_CANCELING, 132 /// The "Minimum Mean Cycle-Canceling" algorithm, which is a 133 /// well-known strongly polynomial method 134 /// \cite goldberg89cyclecanceling. It improves along a 135 /// \ref min_mean_cycle "minimum mean cycle" in each iteration. 136 /// Its running time complexity is \f$O(n^2 m^3 \log n)\f$. 137 MINIMUM_MEAN_CYCLE_CANCELING, 138 /// The "Cancel-and-Tighten" algorithm, which can be viewed as an 139 /// improved version of the previous method 140 /// \cite goldberg89cyclecanceling. 141 /// It is faster both in theory and in practice, its running time 142 /// complexity is \f$O(n^2 m^2 \log n)\f$. 143 CANCEL_AND_TIGHTEN 144 }; 145 146 private: 147 148 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 149 150 typedef std::vector<int> IntVector; 151 typedef std::vector<double> DoubleVector; 152 typedef std::vector<Value> ValueVector; 153 typedef std::vector<Cost> CostVector; 154 typedef std::vector<char> BoolVector; 155 // Note: vector<char> is used instead of vector<bool> for efficiency reasons 156 157 private: 158 159 template <typename KT, typename VT> 160 class StaticVectorMap { 161 public: 162 typedef KT Key; 163 typedef VT Value; 164 StaticVectorMap(std::vector<Value> & v)165 StaticVectorMap(std::vector<Value>& v) : _v(v) {} 166 167 const Value& operator[](const Key& key) const { 168 return _v[StaticDigraph::id(key)]; 169 } 170 171 Value& operator[](const Key& key) { 172 return _v[StaticDigraph::id(key)]; 173 } 174 set(const Key & key,const Value & val)175 void set(const Key& key, const Value& val) { 176 _v[StaticDigraph::id(key)] = val; 177 } 178 179 private: 180 std::vector<Value>& _v; 181 }; 182 183 typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap; 184 typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap; 185 186 private: 187 188 189 // Data related to the underlying digraph 190 const GR &_graph; 191 int _node_num; 192 int _arc_num; 193 int _res_node_num; 194 int _res_arc_num; 195 int _root; 196 197 // Parameters of the problem 198 bool _have_lower; 199 Value _sum_supply; 200 201 // Data structures for storing the digraph 202 IntNodeMap _node_id; 203 IntArcMap _arc_idf; 204 IntArcMap _arc_idb; 205 IntVector _first_out; 206 BoolVector _forward; 207 IntVector _source; 208 IntVector _target; 209 IntVector _reverse; 210 211 // Node and arc data 212 ValueVector _lower; 213 ValueVector _upper; 214 CostVector _cost; 215 ValueVector _supply; 216 217 ValueVector _res_cap; 218 CostVector _pi; 219 220 // Data for a StaticDigraph structure 221 typedef std::pair<int, int> IntPair; 222 StaticDigraph _sgr; 223 std::vector<IntPair> _arc_vec; 224 std::vector<Cost> _cost_vec; 225 IntVector _id_vec; 226 CostArcMap _cost_map; 227 CostNodeMap _pi_map; 228 229 public: 230 231 /// \brief Constant for infinite upper bounds (capacities). 232 /// 233 /// Constant for infinite upper bounds (capacities). 234 /// It is \c std::numeric_limits<Value>::infinity() if available, 235 /// \c std::numeric_limits<Value>::max() otherwise. 236 const Value INF; 237 238 public: 239 240 /// \brief Constructor. 241 /// 242 /// The constructor of the class. 243 /// 244 /// \param graph The digraph the algorithm runs on. CycleCanceling(const GR & graph)245 CycleCanceling(const GR& graph) : 246 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), 247 _cost_map(_cost_vec), _pi_map(_pi), 248 INF(std::numeric_limits<Value>::has_infinity ? 249 std::numeric_limits<Value>::infinity() : 250 std::numeric_limits<Value>::max()) 251 { 252 // Check the number types 253 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, 254 "The flow type of CycleCanceling must be signed"); 255 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 256 "The cost type of CycleCanceling must be signed"); 257 258 // Reset data structures 259 reset(); 260 } 261 262 /// \name Parameters 263 /// The parameters of the algorithm can be specified using these 264 /// functions. 265 266 /// @{ 267 268 /// \brief Set the lower bounds on the arcs. 269 /// 270 /// This function sets the lower bounds on the arcs. 271 /// If it is not used before calling \ref run(), the lower bounds 272 /// will be set to zero on all arcs. 273 /// 274 /// \param map An arc map storing the lower bounds. 275 /// Its \c Value type must be convertible to the \c Value type 276 /// of the algorithm. 277 /// 278 /// \return <tt>(*this)</tt> 279 template <typename LowerMap> lowerMap(const LowerMap & map)280 CycleCanceling& lowerMap(const LowerMap& map) { 281 _have_lower = true; 282 for (ArcIt a(_graph); a != INVALID; ++a) { 283 _lower[_arc_idf[a]] = map[a]; 284 _lower[_arc_idb[a]] = map[a]; 285 } 286 return *this; 287 } 288 289 /// \brief Set the upper bounds (capacities) on the arcs. 290 /// 291 /// This function sets the upper bounds (capacities) on the arcs. 292 /// If it is not used before calling \ref run(), the upper bounds 293 /// will be set to \ref INF on all arcs (i.e. the flow value will be 294 /// unbounded from above). 295 /// 296 /// \param map An arc map storing the upper bounds. 297 /// Its \c Value type must be convertible to the \c Value type 298 /// of the algorithm. 299 /// 300 /// \return <tt>(*this)</tt> 301 template<typename UpperMap> upperMap(const UpperMap & map)302 CycleCanceling& upperMap(const UpperMap& map) { 303 for (ArcIt a(_graph); a != INVALID; ++a) { 304 _upper[_arc_idf[a]] = map[a]; 305 } 306 return *this; 307 } 308 309 /// \brief Set the costs of the arcs. 310 /// 311 /// This function sets the costs of the arcs. 312 /// If it is not used before calling \ref run(), the costs 313 /// will be set to \c 1 on all arcs. 314 /// 315 /// \param map An arc map storing the costs. 316 /// Its \c Value type must be convertible to the \c Cost type 317 /// of the algorithm. 318 /// 319 /// \return <tt>(*this)</tt> 320 template<typename CostMap> costMap(const CostMap & map)321 CycleCanceling& costMap(const CostMap& map) { 322 for (ArcIt a(_graph); a != INVALID; ++a) { 323 _cost[_arc_idf[a]] = map[a]; 324 _cost[_arc_idb[a]] = -map[a]; 325 } 326 return *this; 327 } 328 329 /// \brief Set the supply values of the nodes. 330 /// 331 /// This function sets the supply values of the nodes. 332 /// If neither this function nor \ref stSupply() is used before 333 /// calling \ref run(), the supply of each node will be set to zero. 334 /// 335 /// \param map A node map storing the supply values. 336 /// Its \c Value type must be convertible to the \c Value type 337 /// of the algorithm. 338 /// 339 /// \return <tt>(*this)</tt> 340 template<typename SupplyMap> supplyMap(const SupplyMap & map)341 CycleCanceling& supplyMap(const SupplyMap& map) { 342 for (NodeIt n(_graph); n != INVALID; ++n) { 343 _supply[_node_id[n]] = map[n]; 344 } 345 return *this; 346 } 347 348 /// \brief Set single source and target nodes and a supply value. 349 /// 350 /// This function sets a single source node and a single target node 351 /// and the required flow value. 352 /// If neither this function nor \ref supplyMap() is used before 353 /// calling \ref run(), the supply of each node will be set to zero. 354 /// 355 /// Using this function has the same effect as using \ref supplyMap() 356 /// with a map in which \c k is assigned to \c s, \c -k is 357 /// assigned to \c t and all other nodes have zero supply value. 358 /// 359 /// \param s The source node. 360 /// \param t The target node. 361 /// \param k The required amount of flow from node \c s to node \c t 362 /// (i.e. the supply of \c s and the demand of \c t). 363 /// 364 /// \return <tt>(*this)</tt> stSupply(const Node & s,const Node & t,Value k)365 CycleCanceling& stSupply(const Node& s, const Node& t, Value k) { 366 for (int i = 0; i != _res_node_num; ++i) { 367 _supply[i] = 0; 368 } 369 _supply[_node_id[s]] = k; 370 _supply[_node_id[t]] = -k; 371 return *this; 372 } 373 374 /// @} 375 376 /// \name Execution control 377 /// The algorithm can be executed using \ref run(). 378 379 /// @{ 380 381 /// \brief Run the algorithm. 382 /// 383 /// This function runs the algorithm. 384 /// The paramters can be specified using functions \ref lowerMap(), 385 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). 386 /// For example, 387 /// \code 388 /// CycleCanceling<ListDigraph> cc(graph); 389 /// cc.lowerMap(lower).upperMap(upper).costMap(cost) 390 /// .supplyMap(sup).run(); 391 /// \endcode 392 /// 393 /// This function can be called more than once. All the given parameters 394 /// are kept for the next call, unless \ref resetParams() or \ref reset() 395 /// is used, thus only the modified parameters have to be set again. 396 /// If the underlying digraph was also modified after the construction 397 /// of the class (or the last \ref reset() call), then the \ref reset() 398 /// function must be called. 399 /// 400 /// \param method The cycle-canceling method that will be used. 401 /// For more information, see \ref Method. 402 /// 403 /// \return \c INFEASIBLE if no feasible flow exists, 404 /// \n \c OPTIMAL if the problem has optimal solution 405 /// (i.e. it is feasible and bounded), and the algorithm has found 406 /// optimal flow and node potentials (primal and dual solutions), 407 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost 408 /// and infinite upper bound. It means that the objective function 409 /// is unbounded on that arc, however, note that it could actually be 410 /// bounded over the feasible flows, but this algroithm cannot handle 411 /// these cases. 412 /// 413 /// \see ProblemType, Method 414 /// \see resetParams(), reset() 415 ProblemType run(Method method = CANCEL_AND_TIGHTEN) { 416 ProblemType pt = init(); 417 if (pt != OPTIMAL) return pt; 418 start(method); 419 return OPTIMAL; 420 } 421 422 /// \brief Reset all the parameters that have been given before. 423 /// 424 /// This function resets all the paramaters that have been given 425 /// before using functions \ref lowerMap(), \ref upperMap(), 426 /// \ref costMap(), \ref supplyMap(), \ref stSupply(). 427 /// 428 /// It is useful for multiple \ref run() calls. Basically, all the given 429 /// parameters are kept for the next \ref run() call, unless 430 /// \ref resetParams() or \ref reset() is used. 431 /// If the underlying digraph was also modified after the construction 432 /// of the class or the last \ref reset() call, then the \ref reset() 433 /// function must be used, otherwise \ref resetParams() is sufficient. 434 /// 435 /// For example, 436 /// \code 437 /// CycleCanceling<ListDigraph> cs(graph); 438 /// 439 /// // First run 440 /// cc.lowerMap(lower).upperMap(upper).costMap(cost) 441 /// .supplyMap(sup).run(); 442 /// 443 /// // Run again with modified cost map (resetParams() is not called, 444 /// // so only the cost map have to be set again) 445 /// cost[e] += 100; 446 /// cc.costMap(cost).run(); 447 /// 448 /// // Run again from scratch using resetParams() 449 /// // (the lower bounds will be set to zero on all arcs) 450 /// cc.resetParams(); 451 /// cc.upperMap(capacity).costMap(cost) 452 /// .supplyMap(sup).run(); 453 /// \endcode 454 /// 455 /// \return <tt>(*this)</tt> 456 /// 457 /// \see reset(), run() resetParams()458 CycleCanceling& resetParams() { 459 for (int i = 0; i != _res_node_num; ++i) { 460 _supply[i] = 0; 461 } 462 int limit = _first_out[_root]; 463 for (int j = 0; j != limit; ++j) { 464 _lower[j] = 0; 465 _upper[j] = INF; 466 _cost[j] = _forward[j] ? 1 : -1; 467 } 468 for (int j = limit; j != _res_arc_num; ++j) { 469 _lower[j] = 0; 470 _upper[j] = INF; 471 _cost[j] = 0; 472 _cost[_reverse[j]] = 0; 473 } 474 _have_lower = false; 475 return *this; 476 } 477 478 /// \brief Reset the internal data structures and all the parameters 479 /// that have been given before. 480 /// 481 /// This function resets the internal data structures and all the 482 /// paramaters that have been given before using functions \ref lowerMap(), 483 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). 484 /// 485 /// It is useful for multiple \ref run() calls. Basically, all the given 486 /// parameters are kept for the next \ref run() call, unless 487 /// \ref resetParams() or \ref reset() is used. 488 /// If the underlying digraph was also modified after the construction 489 /// of the class or the last \ref reset() call, then the \ref reset() 490 /// function must be used, otherwise \ref resetParams() is sufficient. 491 /// 492 /// See \ref resetParams() for examples. 493 /// 494 /// \return <tt>(*this)</tt> 495 /// 496 /// \see resetParams(), run() reset()497 CycleCanceling& reset() { 498 // Resize vectors 499 _node_num = countNodes(_graph); 500 _arc_num = countArcs(_graph); 501 _res_node_num = _node_num + 1; 502 _res_arc_num = 2 * (_arc_num + _node_num); 503 _root = _node_num; 504 505 _first_out.resize(_res_node_num + 1); 506 _forward.resize(_res_arc_num); 507 _source.resize(_res_arc_num); 508 _target.resize(_res_arc_num); 509 _reverse.resize(_res_arc_num); 510 511 _lower.resize(_res_arc_num); 512 _upper.resize(_res_arc_num); 513 _cost.resize(_res_arc_num); 514 _supply.resize(_res_node_num); 515 516 _res_cap.resize(_res_arc_num); 517 _pi.resize(_res_node_num); 518 519 _arc_vec.reserve(_res_arc_num); 520 _cost_vec.reserve(_res_arc_num); 521 _id_vec.reserve(_res_arc_num); 522 523 // Copy the graph 524 int i = 0, j = 0, k = 2 * _arc_num + _node_num; 525 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 526 _node_id[n] = i; 527 } 528 i = 0; 529 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 530 _first_out[i] = j; 531 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { 532 _arc_idf[a] = j; 533 _forward[j] = true; 534 _source[j] = i; 535 _target[j] = _node_id[_graph.runningNode(a)]; 536 } 537 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { 538 _arc_idb[a] = j; 539 _forward[j] = false; 540 _source[j] = i; 541 _target[j] = _node_id[_graph.runningNode(a)]; 542 } 543 _forward[j] = false; 544 _source[j] = i; 545 _target[j] = _root; 546 _reverse[j] = k; 547 _forward[k] = true; 548 _source[k] = _root; 549 _target[k] = i; 550 _reverse[k] = j; 551 ++j; ++k; 552 } 553 _first_out[i] = j; 554 _first_out[_res_node_num] = k; 555 for (ArcIt a(_graph); a != INVALID; ++a) { 556 int fi = _arc_idf[a]; 557 int bi = _arc_idb[a]; 558 _reverse[fi] = bi; 559 _reverse[bi] = fi; 560 } 561 562 // Reset parameters 563 resetParams(); 564 return *this; 565 } 566 567 /// @} 568 569 /// \name Query Functions 570 /// The results of the algorithm can be obtained using these 571 /// functions.\n 572 /// The \ref run() function must be called before using them. 573 574 /// @{ 575 576 /// \brief Return the total cost of the found flow. 577 /// 578 /// This function returns the total cost of the found flow. 579 /// Its complexity is O(m). 580 /// 581 /// \note The return type of the function can be specified as a 582 /// template parameter. For example, 583 /// \code 584 /// cc.totalCost<double>(); 585 /// \endcode 586 /// It is useful if the total cost cannot be stored in the \c Cost 587 /// type of the algorithm, which is the default return type of the 588 /// function. 589 /// 590 /// \pre \ref run() must be called before using this function. 591 template <typename Number> totalCost()592 Number totalCost() const { 593 Number c = 0; 594 for (ArcIt a(_graph); a != INVALID; ++a) { 595 int i = _arc_idb[a]; 596 c += static_cast<Number>(_res_cap[i]) * 597 (-static_cast<Number>(_cost[i])); 598 } 599 return c; 600 } 601 602 #ifndef DOXYGEN totalCost()603 Cost totalCost() const { 604 return totalCost<Cost>(); 605 } 606 #endif 607 608 /// \brief Return the flow on the given arc. 609 /// 610 /// This function returns the flow on the given arc. 611 /// 612 /// \pre \ref run() must be called before using this function. flow(const Arc & a)613 Value flow(const Arc& a) const { 614 return _res_cap[_arc_idb[a]]; 615 } 616 617 /// \brief Copy the flow values (the primal solution) into the 618 /// given map. 619 /// 620 /// This function copies the flow value on each arc into the given 621 /// map. The \c Value type of the algorithm must be convertible to 622 /// the \c Value type of the map. 623 /// 624 /// \pre \ref run() must be called before using this function. 625 template <typename FlowMap> flowMap(FlowMap & map)626 void flowMap(FlowMap &map) const { 627 for (ArcIt a(_graph); a != INVALID; ++a) { 628 map.set(a, _res_cap[_arc_idb[a]]); 629 } 630 } 631 632 /// \brief Return the potential (dual value) of the given node. 633 /// 634 /// This function returns the potential (dual value) of the 635 /// given node. 636 /// 637 /// \pre \ref run() must be called before using this function. potential(const Node & n)638 Cost potential(const Node& n) const { 639 return static_cast<Cost>(_pi[_node_id[n]]); 640 } 641 642 /// \brief Copy the potential values (the dual solution) into the 643 /// given map. 644 /// 645 /// This function copies the potential (dual value) of each node 646 /// into the given map. 647 /// The \c Cost type of the algorithm must be convertible to the 648 /// \c Value type of the map. 649 /// 650 /// \pre \ref run() must be called before using this function. 651 template <typename PotentialMap> potentialMap(PotentialMap & map)652 void potentialMap(PotentialMap &map) const { 653 for (NodeIt n(_graph); n != INVALID; ++n) { 654 map.set(n, static_cast<Cost>(_pi[_node_id[n]])); 655 } 656 } 657 658 /// @} 659 660 private: 661 662 // Initialize the algorithm init()663 ProblemType init() { 664 if (_res_node_num <= 1) return INFEASIBLE; 665 666 // Check the sum of supply values 667 _sum_supply = 0; 668 for (int i = 0; i != _root; ++i) { 669 _sum_supply += _supply[i]; 670 } 671 if (_sum_supply > 0) return INFEASIBLE; 672 673 // Check lower and upper bounds 674 LEMON_DEBUG(checkBoundMaps(), 675 "Upper bounds must be greater or equal to the lower bounds"); 676 677 678 // Initialize vectors 679 for (int i = 0; i != _res_node_num; ++i) { 680 _pi[i] = 0; 681 } 682 ValueVector excess(_supply); 683 684 // Remove infinite upper bounds and check negative arcs 685 const Value MAX = std::numeric_limits<Value>::max(); 686 int last_out; 687 if (_have_lower) { 688 for (int i = 0; i != _root; ++i) { 689 last_out = _first_out[i+1]; 690 for (int j = _first_out[i]; j != last_out; ++j) { 691 if (_forward[j]) { 692 Value c = _cost[j] < 0 ? _upper[j] : _lower[j]; 693 if (c >= MAX) return UNBOUNDED; 694 excess[i] -= c; 695 excess[_target[j]] += c; 696 } 697 } 698 } 699 } else { 700 for (int i = 0; i != _root; ++i) { 701 last_out = _first_out[i+1]; 702 for (int j = _first_out[i]; j != last_out; ++j) { 703 if (_forward[j] && _cost[j] < 0) { 704 Value c = _upper[j]; 705 if (c >= MAX) return UNBOUNDED; 706 excess[i] -= c; 707 excess[_target[j]] += c; 708 } 709 } 710 } 711 } 712 Value ex, max_cap = 0; 713 for (int i = 0; i != _res_node_num; ++i) { 714 ex = excess[i]; 715 if (ex < 0) max_cap -= ex; 716 } 717 for (int j = 0; j != _res_arc_num; ++j) { 718 if (_upper[j] >= MAX) _upper[j] = max_cap; 719 } 720 721 // Initialize maps for Circulation and remove non-zero lower bounds 722 ConstMap<Arc, Value> low(0); 723 typedef typename Digraph::template ArcMap<Value> ValueArcMap; 724 typedef typename Digraph::template NodeMap<Value> ValueNodeMap; 725 ValueArcMap cap(_graph), flow(_graph); 726 ValueNodeMap sup(_graph); 727 for (NodeIt n(_graph); n != INVALID; ++n) { 728 sup[n] = _supply[_node_id[n]]; 729 } 730 if (_have_lower) { 731 for (ArcIt a(_graph); a != INVALID; ++a) { 732 int j = _arc_idf[a]; 733 Value c = _lower[j]; 734 cap[a] = _upper[j] - c; 735 sup[_graph.source(a)] -= c; 736 sup[_graph.target(a)] += c; 737 } 738 } else { 739 for (ArcIt a(_graph); a != INVALID; ++a) { 740 cap[a] = _upper[_arc_idf[a]]; 741 } 742 } 743 744 // Find a feasible flow using Circulation 745 Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap> 746 circ(_graph, low, cap, sup); 747 if (!circ.flowMap(flow).run()) return INFEASIBLE; 748 749 // Set residual capacities and handle GEQ supply type 750 if (_sum_supply < 0) { 751 for (ArcIt a(_graph); a != INVALID; ++a) { 752 Value fa = flow[a]; 753 _res_cap[_arc_idf[a]] = cap[a] - fa; 754 _res_cap[_arc_idb[a]] = fa; 755 sup[_graph.source(a)] -= fa; 756 sup[_graph.target(a)] += fa; 757 } 758 for (NodeIt n(_graph); n != INVALID; ++n) { 759 excess[_node_id[n]] = sup[n]; 760 } 761 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { 762 int u = _target[a]; 763 int ra = _reverse[a]; 764 _res_cap[a] = -_sum_supply + 1; 765 _res_cap[ra] = -excess[u]; 766 _cost[a] = 0; 767 _cost[ra] = 0; 768 } 769 } else { 770 for (ArcIt a(_graph); a != INVALID; ++a) { 771 Value fa = flow[a]; 772 _res_cap[_arc_idf[a]] = cap[a] - fa; 773 _res_cap[_arc_idb[a]] = fa; 774 } 775 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { 776 int ra = _reverse[a]; 777 _res_cap[a] = 1; 778 _res_cap[ra] = 0; 779 _cost[a] = 0; 780 _cost[ra] = 0; 781 } 782 } 783 784 return OPTIMAL; 785 } 786 787 // Check if the upper bound is greater or equal to the lower bound 788 // on each arc. checkBoundMaps()789 bool checkBoundMaps() { 790 for (int j = 0; j != _res_arc_num; ++j) { 791 if (_upper[j] < _lower[j]) return false; 792 } 793 return true; 794 } 795 796 // Build a StaticDigraph structure containing the current 797 // residual network buildResidualNetwork()798 void buildResidualNetwork() { 799 _arc_vec.clear(); 800 _cost_vec.clear(); 801 _id_vec.clear(); 802 for (int j = 0; j != _res_arc_num; ++j) { 803 if (_res_cap[j] > 0) { 804 _arc_vec.push_back(IntPair(_source[j], _target[j])); 805 _cost_vec.push_back(_cost[j]); 806 _id_vec.push_back(j); 807 } 808 } 809 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 810 } 811 812 // Execute the algorithm and transform the results start(Method method)813 void start(Method method) { 814 // Execute the algorithm 815 switch (method) { 816 case SIMPLE_CYCLE_CANCELING: 817 startSimpleCycleCanceling(); 818 break; 819 case MINIMUM_MEAN_CYCLE_CANCELING: 820 startMinMeanCycleCanceling(); 821 break; 822 case CANCEL_AND_TIGHTEN: 823 startCancelAndTighten(); 824 break; 825 } 826 827 // Compute node potentials 828 if (method != SIMPLE_CYCLE_CANCELING) { 829 buildResidualNetwork(); 830 typename BellmanFord<StaticDigraph, CostArcMap> 831 ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map); 832 bf.distMap(_pi_map); 833 bf.init(0); 834 bf.start(); 835 } 836 837 // Handle non-zero lower bounds 838 if (_have_lower) { 839 int limit = _first_out[_root]; 840 for (int j = 0; j != limit; ++j) { 841 if (!_forward[j]) _res_cap[j] += _lower[j]; 842 } 843 } 844 } 845 846 // Execute the "Simple Cycle Canceling" method startSimpleCycleCanceling()847 void startSimpleCycleCanceling() { 848 // Constants for computing the iteration limits 849 const int BF_FIRST_LIMIT = 2; 850 const double BF_LIMIT_FACTOR = 1.5; 851 852 typedef StaticVectorMap<StaticDigraph::Arc, Value> FilterMap; 853 typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph; 854 typedef StaticVectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap; 855 typedef typename BellmanFord<ResDigraph, CostArcMap> 856 ::template SetDistMap<CostNodeMap> 857 ::template SetPredMap<PredMap>::Create BF; 858 859 // Build the residual network 860 _arc_vec.clear(); 861 _cost_vec.clear(); 862 for (int j = 0; j != _res_arc_num; ++j) { 863 _arc_vec.push_back(IntPair(_source[j], _target[j])); 864 _cost_vec.push_back(_cost[j]); 865 } 866 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 867 868 FilterMap filter_map(_res_cap); 869 ResDigraph rgr(_sgr, filter_map); 870 std::vector<int> cycle; 871 std::vector<StaticDigraph::Arc> pred(_res_arc_num); 872 PredMap pred_map(pred); 873 BF bf(rgr, _cost_map); 874 bf.distMap(_pi_map).predMap(pred_map); 875 876 int length_bound = BF_FIRST_LIMIT; 877 bool optimal = false; 878 while (!optimal) { 879 bf.init(0); 880 int iter_num = 0; 881 bool cycle_found = false; 882 while (!cycle_found) { 883 // Perform some iterations of the Bellman-Ford algorithm 884 int curr_iter_num = iter_num + length_bound <= _node_num ? 885 length_bound : _node_num - iter_num; 886 iter_num += curr_iter_num; 887 int real_iter_num = curr_iter_num; 888 for (int i = 0; i < curr_iter_num; ++i) { 889 if (bf.processNextWeakRound()) { 890 real_iter_num = i; 891 break; 892 } 893 } 894 if (real_iter_num < curr_iter_num) { 895 // Optimal flow is found 896 optimal = true; 897 break; 898 } else { 899 // Search for node disjoint negative cycles 900 std::vector<int> state(_res_node_num, 0); 901 int id = 0; 902 for (int u = 0; u != _res_node_num; ++u) { 903 if (state[u] != 0) continue; 904 ++id; 905 int v = u; 906 for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ? 907 -1 : rgr.id(rgr.source(pred[v]))) { 908 state[v] = id; 909 } 910 if (v != -1 && state[v] == id) { 911 // A negative cycle is found 912 cycle_found = true; 913 cycle.clear(); 914 StaticDigraph::Arc a = pred[v]; 915 Value d, delta = _res_cap[rgr.id(a)]; 916 cycle.push_back(rgr.id(a)); 917 while (rgr.id(rgr.source(a)) != v) { 918 a = pred_map[rgr.source(a)]; 919 d = _res_cap[rgr.id(a)]; 920 if (d < delta) delta = d; 921 cycle.push_back(rgr.id(a)); 922 } 923 924 // Augment along the cycle 925 for (int i = 0; i < int(cycle.size()); ++i) { 926 int j = cycle[i]; 927 _res_cap[j] -= delta; 928 _res_cap[_reverse[j]] += delta; 929 } 930 } 931 } 932 } 933 934 // Increase iteration limit if no cycle is found 935 if (!cycle_found) { 936 length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR); 937 } 938 } 939 } 940 } 941 942 // Execute the "Minimum Mean Cycle Canceling" method startMinMeanCycleCanceling()943 void startMinMeanCycleCanceling() { 944 typedef Path<StaticDigraph> SPath; 945 typedef typename SPath::ArcIt SPathArcIt; 946 typedef typename HowardMmc<StaticDigraph, CostArcMap> 947 ::template SetPath<SPath>::Create HwMmc; 948 typedef typename HartmannOrlinMmc<StaticDigraph, CostArcMap> 949 ::template SetPath<SPath>::Create HoMmc; 950 951 const double HW_ITER_LIMIT_FACTOR = 1.0; 952 const int HW_ITER_LIMIT_MIN_VALUE = 5; 953 954 const int hw_iter_limit = 955 std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num), 956 HW_ITER_LIMIT_MIN_VALUE); 957 958 SPath cycle; 959 HwMmc hw_mmc(_sgr, _cost_map); 960 hw_mmc.cycle(cycle); 961 buildResidualNetwork(); 962 while (true) { 963 964 typename HwMmc::TerminationCause hw_tc = 965 hw_mmc.findCycleMean(hw_iter_limit); 966 if (hw_tc == HwMmc::ITERATION_LIMIT) { 967 // Howard's algorithm reached the iteration limit, start a 968 // strongly polynomial algorithm instead 969 HoMmc ho_mmc(_sgr, _cost_map); 970 ho_mmc.cycle(cycle); 971 // Find a minimum mean cycle (Hartmann-Orlin algorithm) 972 if (!(ho_mmc.findCycleMean() && ho_mmc.cycleCost() < 0)) break; 973 ho_mmc.findCycle(); 974 } else { 975 // Find a minimum mean cycle (Howard algorithm) 976 if (!(hw_tc == HwMmc::OPTIMAL && hw_mmc.cycleCost() < 0)) break; 977 hw_mmc.findCycle(); 978 } 979 980 // Compute delta value 981 Value delta = INF; 982 for (SPathArcIt a(cycle); a != INVALID; ++a) { 983 Value d = _res_cap[_id_vec[_sgr.id(a)]]; 984 if (d < delta) delta = d; 985 } 986 987 // Augment along the cycle 988 for (SPathArcIt a(cycle); a != INVALID; ++a) { 989 int j = _id_vec[_sgr.id(a)]; 990 _res_cap[j] -= delta; 991 _res_cap[_reverse[j]] += delta; 992 } 993 994 // Rebuild the residual network 995 buildResidualNetwork(); 996 } 997 } 998 999 // Execute the "Cancel-and-Tighten" method startCancelAndTighten()1000 void startCancelAndTighten() { 1001 // Constants for the min mean cycle computations 1002 const double LIMIT_FACTOR = 1.0; 1003 const int MIN_LIMIT = 5; 1004 const double HW_ITER_LIMIT_FACTOR = 1.0; 1005 const int HW_ITER_LIMIT_MIN_VALUE = 5; 1006 1007 const int hw_iter_limit = 1008 std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num), 1009 HW_ITER_LIMIT_MIN_VALUE); 1010 1011 // Contruct auxiliary data vectors 1012 DoubleVector pi(_res_node_num, 0.0); 1013 IntVector level(_res_node_num); 1014 BoolVector reached(_res_node_num); 1015 BoolVector processed(_res_node_num); 1016 IntVector pred_node(_res_node_num); 1017 IntVector pred_arc(_res_node_num); 1018 std::vector<int> stack(_res_node_num); 1019 std::vector<int> proc_vector(_res_node_num); 1020 1021 // Initialize epsilon 1022 double epsilon = 0; 1023 for (int a = 0; a != _res_arc_num; ++a) { 1024 if (_res_cap[a] > 0 && -_cost[a] > epsilon) 1025 epsilon = -_cost[a]; 1026 } 1027 1028 // Start phases 1029 Tolerance<double> tol; 1030 tol.epsilon(1e-6); 1031 int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num))); 1032 if (limit < MIN_LIMIT) limit = MIN_LIMIT; 1033 int iter = limit; 1034 while (epsilon * _res_node_num >= 1) { 1035 // Find and cancel cycles in the admissible network using DFS 1036 for (int u = 0; u != _res_node_num; ++u) { 1037 reached[u] = false; 1038 processed[u] = false; 1039 } 1040 int stack_head = -1; 1041 int proc_head = -1; 1042 for (int start = 0; start != _res_node_num; ++start) { 1043 if (reached[start]) continue; 1044 1045 // New start node 1046 reached[start] = true; 1047 pred_arc[start] = -1; 1048 pred_node[start] = -1; 1049 1050 // Find the first admissible outgoing arc 1051 double p = pi[start]; 1052 int a = _first_out[start]; 1053 int last_out = _first_out[start+1]; 1054 for (; a != last_out && (_res_cap[a] == 0 || 1055 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; 1056 if (a == last_out) { 1057 processed[start] = true; 1058 proc_vector[++proc_head] = start; 1059 continue; 1060 } 1061 stack[++stack_head] = a; 1062 1063 while (stack_head >= 0) { 1064 int sa = stack[stack_head]; 1065 int u = _source[sa]; 1066 int v = _target[sa]; 1067 1068 if (!reached[v]) { 1069 // A new node is reached 1070 reached[v] = true; 1071 pred_node[v] = u; 1072 pred_arc[v] = sa; 1073 p = pi[v]; 1074 a = _first_out[v]; 1075 last_out = _first_out[v+1]; 1076 for (; a != last_out && (_res_cap[a] == 0 || 1077 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; 1078 stack[++stack_head] = a == last_out ? -1 : a; 1079 } else { 1080 if (!processed[v]) { 1081 // A cycle is found 1082 int n, w = u; 1083 Value d, delta = _res_cap[sa]; 1084 for (n = u; n != v; n = pred_node[n]) { 1085 d = _res_cap[pred_arc[n]]; 1086 if (d <= delta) { 1087 delta = d; 1088 w = pred_node[n]; 1089 } 1090 } 1091 1092 // Augment along the cycle 1093 _res_cap[sa] -= delta; 1094 _res_cap[_reverse[sa]] += delta; 1095 for (n = u; n != v; n = pred_node[n]) { 1096 int pa = pred_arc[n]; 1097 _res_cap[pa] -= delta; 1098 _res_cap[_reverse[pa]] += delta; 1099 } 1100 for (n = u; stack_head > 0 && n != w; n = pred_node[n]) { 1101 --stack_head; 1102 reached[n] = false; 1103 } 1104 u = w; 1105 } 1106 v = u; 1107 1108 // Find the next admissible outgoing arc 1109 p = pi[v]; 1110 a = stack[stack_head] + 1; 1111 last_out = _first_out[v+1]; 1112 for (; a != last_out && (_res_cap[a] == 0 || 1113 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; 1114 stack[stack_head] = a == last_out ? -1 : a; 1115 } 1116 1117 while (stack_head >= 0 && stack[stack_head] == -1) { 1118 processed[v] = true; 1119 proc_vector[++proc_head] = v; 1120 if (--stack_head >= 0) { 1121 // Find the next admissible outgoing arc 1122 v = _source[stack[stack_head]]; 1123 p = pi[v]; 1124 a = stack[stack_head] + 1; 1125 last_out = _first_out[v+1]; 1126 for (; a != last_out && (_res_cap[a] == 0 || 1127 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; 1128 stack[stack_head] = a == last_out ? -1 : a; 1129 } 1130 } 1131 } 1132 } 1133 1134 // Tighten potentials and epsilon 1135 if (--iter > 0) { 1136 for (int u = 0; u != _res_node_num; ++u) { 1137 level[u] = 0; 1138 } 1139 for (int i = proc_head; i > 0; --i) { 1140 int u = proc_vector[i]; 1141 double p = pi[u]; 1142 int l = level[u] + 1; 1143 int last_out = _first_out[u+1]; 1144 for (int a = _first_out[u]; a != last_out; ++a) { 1145 int v = _target[a]; 1146 if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) && 1147 l > level[v]) level[v] = l; 1148 } 1149 } 1150 1151 // Modify potentials 1152 double q = std::numeric_limits<double>::max(); 1153 for (int u = 0; u != _res_node_num; ++u) { 1154 int lu = level[u]; 1155 double p, pu = pi[u]; 1156 int last_out = _first_out[u+1]; 1157 for (int a = _first_out[u]; a != last_out; ++a) { 1158 if (_res_cap[a] == 0) continue; 1159 int v = _target[a]; 1160 int ld = lu - level[v]; 1161 if (ld > 0) { 1162 p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1); 1163 if (p < q) q = p; 1164 } 1165 } 1166 } 1167 for (int u = 0; u != _res_node_num; ++u) { 1168 pi[u] -= q * level[u]; 1169 } 1170 1171 // Modify epsilon 1172 epsilon = 0; 1173 for (int u = 0; u != _res_node_num; ++u) { 1174 double curr, pu = pi[u]; 1175 int last_out = _first_out[u+1]; 1176 for (int a = _first_out[u]; a != last_out; ++a) { 1177 if (_res_cap[a] == 0) continue; 1178 curr = _cost[a] + pu - pi[_target[a]]; 1179 if (-curr > epsilon) epsilon = -curr; 1180 } 1181 } 1182 } else { 1183 typedef HowardMmc<StaticDigraph, CostArcMap> HwMmc; 1184 typedef HartmannOrlinMmc<StaticDigraph, CostArcMap> HoMmc; 1185 typedef typename BellmanFord<StaticDigraph, CostArcMap> 1186 ::template SetDistMap<CostNodeMap>::Create BF; 1187 1188 // Set epsilon to the minimum cycle mean 1189 Cost cycle_cost = 0; 1190 int cycle_size = 1; 1191 buildResidualNetwork(); 1192 HwMmc hw_mmc(_sgr, _cost_map); 1193 if (hw_mmc.findCycleMean(hw_iter_limit) == HwMmc::ITERATION_LIMIT) { 1194 // Howard's algorithm reached the iteration limit, start a 1195 // strongly polynomial algorithm instead 1196 HoMmc ho_mmc(_sgr, _cost_map); 1197 ho_mmc.findCycleMean(); 1198 epsilon = -ho_mmc.cycleMean(); 1199 cycle_cost = ho_mmc.cycleCost(); 1200 cycle_size = ho_mmc.cycleSize(); 1201 } else { 1202 // Set epsilon 1203 epsilon = -hw_mmc.cycleMean(); 1204 cycle_cost = hw_mmc.cycleCost(); 1205 cycle_size = hw_mmc.cycleSize(); 1206 } 1207 1208 // Compute feasible potentials for the current epsilon 1209 for (int i = 0; i != int(_cost_vec.size()); ++i) { 1210 _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost; 1211 } 1212 BF bf(_sgr, _cost_map); 1213 bf.distMap(_pi_map); 1214 bf.init(0); 1215 bf.start(); 1216 for (int u = 0; u != _res_node_num; ++u) { 1217 pi[u] = static_cast<double>(_pi[u]) / cycle_size; 1218 } 1219 1220 iter = limit; 1221 } 1222 } 1223 } 1224 1225 }; //class CycleCanceling 1226 1227 ///@} 1228 1229 } //namespace lemon 1230 1231 #endif //LEMON_CYCLE_CANCELING_H 1232