Changeset 2556:74c2c81055e1 in lemon-0.x for lemon/network_simplex.h
- Timestamp:
- 01/13/08 11:32:14 (16 years ago)
- Branch:
- default
- Phase:
- public
- Convert:
- svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3437
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/network_simplex.h
r2553 r2556 23 23 /// 24 24 /// \file 25 /// \brief The network simplex algorithm for finding a minimum cost 26 /// flow. 25 /// \brief The network simplex algorithm for finding a minimum cost flow. 27 26 28 27 #include <limits> … … 41 40 42 41 43 // / \briefState constant for edges at their lower bounds.44 #define LOWER 45 // / \briefState constant for edges in the spanning tree.46 #define TREE 47 // / \briefState constant for edges at their upper bounds.48 #define UPPER 42 // State constant for edges at their lower bounds. 43 #define LOWER 1 44 // State constant for edges in the spanning tree. 45 #define TREE 0 46 // State constant for edges at their upper bounds. 47 #define UPPER -1 49 48 50 49 #ifdef EDGE_BLOCK_PIVOT 51 50 #include <cmath> 52 /// \brief Lower bound for the size of blocks. 53 #define MIN_BLOCK_SIZE 10 51 #define MIN_BLOCK_SIZE 10 54 52 #endif 55 53 56 54 #ifdef CANDIDATE_LIST_PIVOT 57 55 #include <vector> 58 #define LIST_LENGTH_DIV 59 #define MINOR_LIMIT_DIV 56 #define LIST_LENGTH_DIV 50 57 #define MINOR_LIMIT_DIV 200 60 58 #endif 61 59 … … 64 62 #include <algorithm> 65 63 #define LIST_LENGTH_DIV 100 66 #define LOWER_DIV 64 #define LOWER_DIV 4 67 65 #endif 68 66 … … 75 73 /// finding a minimum cost flow. 76 74 /// 77 /// \ref lemon::NetworkSimplex "NetworkSimplex" implements the78 /// network simplex algorithm forfinding a minimum cost flow.75 /// \ref NetworkSimplex implements the network simplex algorithm for 76 /// finding a minimum cost flow. 79 77 /// 80 78 /// \param Graph The directed graph type the algorithm runs on. … … 85 83 /// 86 84 /// \warning 87 /// - Edge capacities and costs should be non negative integers.88 /// 85 /// - Edge capacities and costs should be non-negative integers. 86 /// However \c CostMap::Value should be signed type. 89 87 /// - Supply values should be signed integers. 90 88 /// - \c LowerMap::Value must be convertible to 91 /// 92 /// 89 /// \c CapacityMap::Value and \c CapacityMap::Value must be 90 /// convertible to \c SupplyMap::Value. 93 91 /// 94 92 /// \author Peter Kovacs … … 108 106 109 107 typedef SmartGraph SGraph; 110 typedef typename SGraph::Node Node; 111 typedef typename SGraph::NodeIt NodeIt; 112 typedef typename SGraph::Edge Edge; 113 typedef typename SGraph::EdgeIt EdgeIt; 114 typedef typename SGraph::InEdgeIt InEdgeIt; 115 typedef typename SGraph::OutEdgeIt OutEdgeIt; 108 GRAPH_TYPEDEFS(typename SGraph); 116 109 117 110 typedef typename SGraph::template EdgeMap<Lower> SLowerMap; … … 132 125 public: 133 126 134 /// \briefThe type of the flow map.127 /// The type of the flow map. 135 128 typedef typename Graph::template EdgeMap<Capacity> FlowMap; 136 /// \briefThe type of the potential map.129 /// The type of the potential map. 137 130 typedef typename Graph::template NodeMap<Cost> PotentialMap; 138 131 139 132 protected: 140 133 141 /// \briefMap adaptor class for handling reduced edge costs.134 /// Map adaptor class for handling reduced edge costs. 142 135 class ReducedCostMap : public MapBase<Edge, Cost> 143 136 { … … 151 144 152 145 ReducedCostMap( const SGraph &_gr, 153 154 155 146 const SCostMap &_cm, 147 const SPotentialMap &_pm ) : 148 gr(_gr), cost_map(_cm), pot_map(_pm) {} 156 149 157 150 Cost operator[](const Edge &e) const { 158 151 return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)]; 159 152 } 160 153 … … 163 156 protected: 164 157 165 /// \briefThe directed graph the algorithm runs on.158 /// The directed graph the algorithm runs on. 166 159 SGraph graph; 167 /// \briefThe original graph.160 /// The original graph. 168 161 const Graph &graph_ref; 169 /// \briefThe original lower bound map.162 /// The original lower bound map. 170 163 const LowerMap *lower; 171 /// \briefThe capacity map.164 /// The capacity map. 172 165 SCapacityMap capacity; 173 /// \briefThe cost map.166 /// The cost map. 174 167 SCostMap cost; 175 /// \briefThe supply map.168 /// The supply map. 176 169 SSupplyMap supply; 177 /// \briefThe reduced cost map.170 /// The reduced cost map. 178 171 ReducedCostMap red_cost; 179 /// \brief The sum of supply values equals zero.180 172 bool valid_supply; 181 173 182 /// \briefThe edge map of the current flow.174 /// The edge map of the current flow. 183 175 SCapacityMap flow; 184 /// \brief The edge map of the found flow on the original graph. 176 /// The potential node map. 177 SPotentialMap potential; 185 178 FlowMap flow_result; 186 /// \brief The potential node map.187 SPotentialMap potential;188 /// \brief The potential node map on the original graph.189 179 PotentialMap potential_result; 190 180 191 /// \briefNode reference for the original graph.181 /// Node reference for the original graph. 192 182 NodeRefMap node_ref; 193 /// \briefEdge reference for the original graph.183 /// Edge reference for the original graph. 194 184 EdgeRefMap edge_ref; 195 185 196 /// \brief Thedepth node map of the spanning tree structure.186 /// The \c depth node map of the spanning tree structure. 197 187 IntNodeMap depth; 198 /// \brief Theparent node map of the spanning tree structure.188 /// The \c parent node map of the spanning tree structure. 199 189 NodeNodeMap parent; 200 /// \brief Thepred_edge node map of the spanning tree structure.190 /// The \c pred_edge node map of the spanning tree structure. 201 191 EdgeNodeMap pred_edge; 202 /// \brief Thethread node map of the spanning tree structure.192 /// The \c thread node map of the spanning tree structure. 203 193 NodeNodeMap thread; 204 /// \brief Theforward node map of the spanning tree structure.194 /// The \c forward node map of the spanning tree structure. 205 195 BoolNodeMap forward; 206 /// \briefThe state edge map.196 /// The state edge map. 207 197 IntEdgeMap state; 208 198 /// The root node of the starting spanning tree. 199 Node root; 200 201 // The entering edge of the current pivot iteration. 202 Edge in_edge; 203 // Temporary nodes used in the current pivot iteration. 204 Node join, u_in, v_in, u_out, v_out; 205 Node right, first, second, last; 206 Node stem, par_stem, new_stem; 207 // The maximum augment amount along the found cycle in the current 208 // pivot iteration. 209 Capacity delta; 209 210 210 211 #ifdef EDGE_BLOCK_PIVOT 211 /// \briefThe size of blocks for the "Edge Block" pivot rule.212 /// The size of blocks for the "Edge Block" pivot rule. 212 213 int block_size; 213 214 #endif … … 235 236 #endif 236 237 237 // Root node of the starting spanning tree.238 Node root;239 // The entering edge of the current pivot iteration.240 Edge in_edge;241 // Temporary nodes used in the current pivot iteration.242 Node join, u_in, v_in, u_out, v_out;243 Node right, first, second, last;244 Node stem, par_stem, new_stem;245 // The maximum augment amount along the cycle in the current pivot246 // iteration.247 Capacity delta;248 249 238 public : 250 239 … … 259 248 /// \param _supply The supply values of the nodes (signed). 260 249 NetworkSimplex( const Graph &_graph, 261 262 263 264 250 const LowerMap &_lower, 251 const CapacityMap &_capacity, 252 const CostMap &_cost, 253 const SupplyMap &_supply ) : 265 254 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), 266 255 supply(graph), flow(graph), flow_result(_graph), potential(graph), … … 273 262 Supply sum = 0; 274 263 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) 275 264 sum += _supply[n]; 276 265 if (!(valid_supply = sum == 0)) return; 277 266 … … 280 269 graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref)); 281 270 copyGraph(graph, graph_ref) 282 283 284 285 286 287 // Removing non zero lower bounds271 .edgeMap(cost, _cost) 272 .nodeRef(node_ref) 273 .edgeRef(edge_ref) 274 .run(); 275 276 // Removing non-zero lower bounds 288 277 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { 289 278 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; 290 279 } 291 280 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { 292 293 294 295 296 297 281 Supply s = _supply[n]; 282 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) 283 s += _lower[e]; 284 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) 285 s -= _lower[e]; 286 supply[node_ref[n]] = s; 298 287 } 299 288 } … … 308 297 /// \param _supply The supply values of the nodes (signed). 309 298 NetworkSimplex( const Graph &_graph, 310 311 312 299 const CapacityMap &_capacity, 300 const CostMap &_cost, 301 const SupplyMap &_supply ) : 313 302 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), 314 303 supply(graph), flow(graph), flow_result(_graph), potential(graph), … … 321 310 Supply sum = 0; 322 311 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) 323 312 sum += _supply[n]; 324 313 if (!(valid_supply = sum == 0)) return; 325 314 326 315 // Copying graph_ref to graph 327 316 copyGraph(graph, graph_ref) 328 329 330 331 332 333 317 .edgeMap(capacity, _capacity) 318 .edgeMap(cost, _cost) 319 .nodeMap(supply, _supply) 320 .nodeRef(node_ref) 321 .edgeRef(edge_ref) 322 .run(); 334 323 } 335 324 … … 345 334 /// \param _t The target node. 346 335 /// \param _flow_value The required amount of flow from node \c _s 347 /// to node \c _t (i.e. the supply of \c _s and the demand of 348 /// \c _t). 336 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). 349 337 NetworkSimplex( const Graph &_graph, 350 351 352 353 354 355 338 const LowerMap &_lower, 339 const CapacityMap &_capacity, 340 const CostMap &_cost, 341 typename Graph::Node _s, 342 typename Graph::Node _t, 343 typename SupplyMap::Value _flow_value ) : 356 344 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), 357 345 supply(graph), flow(graph), flow_result(_graph), potential(graph), … … 363 351 // Copying graph_ref to graph 364 352 copyGraph(graph, graph_ref) 365 366 367 368 369 370 // Removing non zero lower bounds353 .edgeMap(cost, _cost) 354 .nodeRef(node_ref) 355 .edgeRef(edge_ref) 356 .run(); 357 358 // Removing non-zero lower bounds 371 359 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { 372 360 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; 373 361 } 374 362 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { 375 376 377 378 379 380 381 382 363 Supply s = 0; 364 if (n == _s) s = _flow_value; 365 if (n == _t) s = -_flow_value; 366 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) 367 s += _lower[e]; 368 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) 369 s -= _lower[e]; 370 supply[node_ref[n]] = s; 383 371 } 384 372 valid_supply = true; … … 395 383 /// \param _t The target node. 396 384 /// \param _flow_value The required amount of flow from node \c _s 397 /// to node \c _t (i.e. the supply of \c _s and the demand of 398 /// \c _t). 385 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). 399 386 NetworkSimplex( const Graph &_graph, 400 401 402 403 404 387 const CapacityMap &_capacity, 388 const CostMap &_cost, 389 typename Graph::Node _s, 390 typename Graph::Node _t, 391 typename SupplyMap::Value _flow_value ) : 405 392 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), 406 393 supply(graph, 0), flow(graph), flow_result(_graph), potential(graph), … … 412 399 // Copying graph_ref to graph 413 400 copyGraph(graph, graph_ref) 414 415 416 417 418 401 .edgeMap(capacity, _capacity) 402 .edgeMap(cost, _cost) 403 .nodeRef(node_ref) 404 .edgeRef(edge_ref) 405 .run(); 419 406 supply[node_ref[_s]] = _flow_value; 420 407 supply[node_ref[_t]] = -_flow_value; 421 408 valid_supply = true; 409 } 410 411 /// \brief Runs the algorithm. 412 /// 413 /// Runs the algorithm. 414 /// 415 /// \return \c true if a feasible flow can be found. 416 bool run() { 417 return init() && start(); 422 418 } 423 419 … … 451 447 Cost c = 0; 452 448 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) 453 449 c += flow_result[e] * cost[edge_ref[e]]; 454 450 return c; 455 }456 457 /// \brief Runs the algorithm.458 ///459 /// Runs the algorithm.460 ///461 /// \return \c true if a feasible flow can be found.462 bool run() {463 return init() && start();464 451 } 465 452 … … 473 460 // Initializing state and flow maps 474 461 for (EdgeIt e(graph); e != INVALID; ++e) { 475 476 462 flow[e] = 0; 463 state[e] = LOWER; 477 464 } 478 465 … … 492 479 Cost max_cost = std::numeric_limits<Cost>::max() / 4; 493 480 for (NodeIt u(graph); u != INVALID; ++u) { 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 481 if (u == root) continue; 482 thread[last] = u; 483 last = u; 484 parent[u] = root; 485 depth[u] = 1; 486 sum += supply[u]; 487 if (supply[u] >= 0) { 488 e = graph.addEdge(u, root); 489 flow[e] = supply[u]; 490 forward[u] = true; 491 potential[u] = max_cost; 492 } else { 493 e = graph.addEdge(root, u); 494 flow[e] = -supply[u]; 495 forward[u] = false; 496 potential[u] = -max_cost; 497 } 498 cost[e] = max_cost; 499 capacity[e] = std::numeric_limits<Capacity>::max(); 500 state[e] = TREE; 501 pred_edge[u] = e; 515 502 } 516 503 thread[last] = root; … … 521 508 block_size = 2 * int(sqrt(countEdges(graph))); 522 509 if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE; 523 // block_size = edge_num >= BLOCK_NUM * MIN_BLOCK_SIZE ?524 // edge_num / BLOCK_NUM : MIN_BLOCK_SIZE;525 510 #endif 526 511 #ifdef CANDIDATE_LIST_PIVOT … … 544 529 bool findEnteringEdge(EdgeIt &next_edge) { 545 530 for (EdgeIt e = next_edge; e != INVALID; ++e) { 546 547 548 549 550 531 if (state[e] * red_cost[e] < 0) { 532 in_edge = e; 533 next_edge = ++e; 534 return true; 535 } 551 536 } 552 537 for (EdgeIt e(graph); e != next_edge; ++e) { 553 554 555 556 557 538 if (state[e] * red_cost[e] < 0) { 539 in_edge = e; 540 next_edge = ++e; 541 return true; 542 } 558 543 } 559 544 return false; … … 567 552 Cost min = 0; 568 553 for (EdgeIt e(graph); e != INVALID; ++e) { 569 570 571 572 554 if (state[e] * red_cost[e] < min) { 555 min = state[e] * red_cost[e]; 556 in_edge = e; 557 } 573 558 } 574 559 return min < 0; … … 585 570 int cnt = 0; 586 571 for (EdgeIt e = next_edge; e != INVALID; ++e) { 587 588 589 590 591 592 593 594 572 if ((curr = state[e] * red_cost[e]) < min) { 573 min = curr; 574 min_edge = e; 575 } 576 if (++cnt == block_size) { 577 if (min < 0) break; 578 cnt = 0; 579 } 595 580 } 596 581 if (!(min < 0)) { 597 598 599 600 601 602 603 604 605 606 582 for (EdgeIt e(graph); e != next_edge; ++e) { 583 if ((curr = state[e] * red_cost[e]) < min) { 584 min = curr; 585 min_edge = e; 586 } 587 if (++cnt == block_size) { 588 if (min < 0) break; 589 cnt = 0; 590 } 591 } 607 592 } 608 593 in_edge = min_edge; 609 594 if ((next_edge = ++min_edge) == INVALID) 610 595 next_edge = EdgeIt(graph); 611 596 return min < 0; 612 597 } … … 620 605 621 606 if (minor_count >= minor_limit || candidates.size() == 0) { 622 623 624 625 626 627 628 629 630 607 // Major iteration 608 candidates.clear(); 609 for (EdgeIt e(graph); e != INVALID; ++e) { 610 if (state[e] * red_cost[e] < 0) { 611 candidates.push_back(e); 612 if (candidates.size() == list_length) break; 613 } 614 } 615 if (candidates.size() == 0) return false; 631 616 } 632 617 … … 637 622 for (int i = 0; i < candidates.size(); ++i) { 638 623 e = candidates[i]; 639 640 641 642 624 if (state[e] * red_cost[e] < min) { 625 min = state[e] * red_cost[e]; 626 in_edge = e; 627 } 643 628 } 644 629 return true; … … 656 641 public: 657 642 SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) : 658 643 st(_st), rc(_rc) {} 659 644 bool operator()(const Edge &e1, const Edge &e2) { 660 645 return st[e1] * rc[e1] < st[e2] * rc[e2]; 661 646 } 662 647 }; … … 669 654 // Minor iteration 670 655 while (list_index < candidates.size()) { 671 672 656 in_edge = candidates[list_index++]; 657 if (state[in_edge] * red_cost[in_edge] < 0) return true; 673 658 } 674 659 … … 677 662 Cost curr, min = 0; 678 663 for (EdgeIt e(graph); e != INVALID; ++e) { 679 680 681 682 683 664 if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) { 665 candidates.push_back(e); 666 if (curr < min) min = curr; 667 if (candidates.size() == list_length) break; 668 } 684 669 } 685 670 if (candidates.size() == 0) return false; … … 697 682 Node v = graph.target(in_edge); 698 683 while (u != v) { 699 700 701 702 703 704 684 if (depth[u] == depth[v]) { 685 u = parent[u]; 686 v = parent[v]; 687 } 688 else if (depth[u] > depth[v]) u = parent[u]; 689 else v = parent[v]; 705 690 } 706 691 return u; … … 713 698 // of the cycle 714 699 if (state[in_edge] == LOWER) { 715 716 second= graph.target(in_edge);700 first = graph.source(in_edge); 701 second = graph.target(in_edge); 717 702 } else { 718 719 second= graph.source(in_edge);703 first = graph.target(in_edge); 704 second = graph.source(in_edge); 720 705 } 721 706 delta = capacity[in_edge]; … … 727 712 // root node 728 713 for (Node u = first; u != join; u = parent[u]) { 729 730 731 732 733 734 735 736 737 714 e = pred_edge[u]; 715 d = forward[u] ? flow[e] : capacity[e] - flow[e]; 716 if (d < delta) { 717 delta = d; 718 u_out = u; 719 u_in = first; 720 v_in = second; 721 result = true; 722 } 738 723 } 739 724 // Searching the cycle along the path form the second node to the 740 725 // root node 741 726 for (Node u = second; u != join; u = parent[u]) { 742 743 744 745 746 747 748 749 750 727 e = pred_edge[u]; 728 d = forward[u] ? capacity[e] - flow[e] : flow[e]; 729 if (d <= delta) { 730 delta = d; 731 u_out = u; 732 u_in = second; 733 v_in = first; 734 result = true; 735 } 751 736 } 752 737 return result; 753 738 } 754 739 755 /// \brief Changes flow andstate edge maps.740 /// \brief Changes \c flow and \c state edge maps. 756 741 void changeFlows(bool change) { 757 742 // Augmenting along the cycle 758 743 if (delta > 0) { 759 760 761 762 763 764 765 766 744 Capacity val = state[in_edge] * delta; 745 flow[in_edge] += val; 746 for (Node u = graph.source(in_edge); u != join; u = parent[u]) { 747 flow[pred_edge[u]] += forward[u] ? -val : val; 748 } 749 for (Node u = graph.target(in_edge); u != join; u = parent[u]) { 750 flow[pred_edge[u]] += forward[u] ? val : -val; 751 } 767 752 } 768 753 // Updating the state of the entering and leaving edges 769 754 if (change) { 770 771 772 755 state[in_edge] = TREE; 756 state[pred_edge[u_out]] = 757 (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER; 773 758 } else { 774 775 } 776 } 777 778 /// \brief Updates thread andparent node maps.759 state[in_edge] = -state[in_edge]; 760 } 761 } 762 763 /// \brief Updates \c thread and \c parent node maps. 779 764 void updateThreadParent() { 780 765 Node u; … … 784 769 bool par_first = false; 785 770 if (join == v_out) { 786 787 788 789 790 791 771 for (u = join; u != u_in && u != v_in; u = thread[u]) ; 772 if (u == v_in) { 773 par_first = true; 774 while (thread[u] != u_out) u = thread[u]; 775 first = u; 776 } 792 777 } 793 778 … … 797 782 right = thread[u]; 798 783 if (thread[v_in] == u_out) { 799 800 784 for (last = u; depth[last] > depth[u_out]; last = thread[last]) ; 785 if (last == u_out) last = thread[last]; 801 786 } 802 787 else last = thread[v_in]; … … 806 791 par_stem = v_in; 807 792 while (stem != u_out) { 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 793 thread[u] = new_stem = parent[stem]; 794 795 // Finding the node just before the stem node (u) according to 796 // the original thread index 797 for (u = new_stem; thread[u] != stem; u = thread[u]) ; 798 thread[u] = right; 799 800 // Changing the parent node of stem and shifting stem and 801 // par_stem nodes 802 parent[stem] = par_stem; 803 par_stem = stem; 804 stem = new_stem; 805 806 // Finding the last successor of stem (u) and the node after it 807 // (right) according to the thread index 808 for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ; 809 right = thread[u]; 825 810 } 826 811 parent[u_out] = par_stem; … … 828 813 829 814 if (join == v_out && par_first) { 830 815 if (first != v_in) thread[first] = right; 831 816 } else { 832 833 834 } 835 } 836 837 /// \brief Updates pred_edge andforward node maps.817 for (u = v_out; thread[u] != u_out; u = thread[u]) ; 818 thread[u] = right; 819 } 820 } 821 822 /// \brief Updates \c pred_edge and \c forward node maps. 838 823 void updatePredEdge() { 839 824 Node u = u_out, v; 840 825 while (u != u_in) { 841 842 843 844 826 v = parent[u]; 827 pred_edge[u] = pred_edge[v]; 828 forward[u] = !forward[v]; 829 u = v; 845 830 } 846 831 pred_edge[u_in] = in_edge; … … 848 833 } 849 834 850 /// \brief Updates depth andpotential node maps.835 /// \brief Updates \c depth and \c potential node maps. 851 836 void updateDepthPotential() { 852 837 depth[u_in] = depth[v_in] + 1; 853 838 potential[u_in] = forward[u_in] ? 854 855 839 potential[v_in] + cost[pred_edge[u_in]] : 840 potential[v_in] - cost[pred_edge[u_in]]; 856 841 857 842 Node u = thread[u_in], v; 858 843 while (true) { 859 860 861 862 863 864 865 866 844 v = parent[u]; 845 if (v == INVALID) break; 846 depth[u] = depth[v] + 1; 847 potential[u] = forward[u] ? 848 potential[v] + cost[pred_edge[u]] : 849 potential[v] - cost[pred_edge[u]]; 850 if (depth[u] <= depth[v_in]) break; 851 u = thread[u]; 867 852 } 868 853 } … … 881 866 #endif 882 867 { 883 884 885 886 887 888 889 890 868 join = findJoinNode(); 869 bool change = findLeavingEdge(); 870 changeFlows(change); 871 if (change) { 872 updateThreadParent(); 873 updatePredEdge(); 874 updateDepthPotential(); 875 } 891 876 #ifdef _DEBUG_ITER_ 892 877 ++iter_num; 893 878 #endif 894 879 } … … 896 881 #ifdef _DEBUG_ITER_ 897 882 std::cout << "Network Simplex algorithm finished. " << iter_num 898 883 << " pivot iterations performed." << std::endl; 899 884 #endif 900 885 … … 902 887 // artificial edges 903 888 for (InEdgeIt e(graph, root); e != INVALID; ++e) 904 889 if (flow[e] > 0) return false; 905 890 for (OutEdgeIt e(graph, root); e != INVALID; ++e) 906 891 if (flow[e] > 0) return false; 907 892 908 893 // Copying flow values to flow_result 909 894 if (lower) { 910 911 895 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) 896 flow_result[e] = (*lower)[e] + flow[edge_ref[e]]; 912 897 } else { 913 914 898 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) 899 flow_result[e] = flow[edge_ref[e]]; 915 900 } 916 901 // Copying potential values to potential_result 917 902 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) 918 903 potential_result[n] = potential[node_ref[n]]; 919 904 920 905 return true;
Note: See TracChangeset
for help on using the changeset viewer.