Changeset 2625:c51b320bc51c in lemon0.x for lemon/cost_scaling.h
 Timestamp:
 10/15/08 14:04:11 (14 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@3510
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/cost_scaling.h
r2623 r2625 21 21 22 22 /// \ingroup min_cost_flow 23 ///24 23 /// \file 25 24 /// \brief Cost scaling algorithm for finding a minimum cost flow. … … 43 42 /// 44 43 /// \ref CostScaling implements the cost scaling algorithm performing 45 /// generalized pushrelabel operations for finding a minimum cost44 /// augment/push and relabel operations for finding a minimum cost 46 45 /// flow. 47 46 /// … … 117 116 118 117 ///\e 119 typename Map::Value operator[](const ResEdge &e) const {120 return ResGraph::forward(e) ? 118 inline typename Map::Value operator[](const ResEdge &e) const { 119 return ResGraph::forward(e) ? _cost_map[e] : _cost_map[e]; 121 120 } 122 121 … … 143 142 144 143 ///\e 145 LCost operator[](const Edge &e) const {144 inline LCost operator[](const Edge &e) const { 146 145 return _cost_map[e] + _pot_map[_gr.source(e)] 147 146  _pot_map[_gr.target(e)]; … … 149 148 150 149 }; //class ReducedCostMap 151 152 private:153 154 // Scaling factor155 static const int ALPHA = 4;156 157 // Paramters for heuristics158 static const int BF_HEURISTIC_EPSILON_BOUND = 5000;159 static const int BF_HEURISTIC_BOUND_FACTOR = 3;160 150 161 151 private: … … 192 182 // The epsilon parameter used for cost scaling 193 183 LCost _epsilon; 184 // The scaling factor 185 int _alpha; 194 186 195 187 public: … … 214 206 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0) 215 207 { 216 // Remov ingnonzero lower bounds208 // Remove nonzero lower bounds 217 209 _capacity = subMap(capacity, lower); 218 210 Supply sum = 0; … … 246 238 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0) 247 239 { 248 // Check ingthe sum of supply values240 // Check the sum of supply values 249 241 Supply sum = 0; 250 242 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; … … 275 267 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0) 276 268 { 277 // Remov ingnonzero lower bounds269 // Remove nonzero lower bounds 278 270 _capacity = subMap(capacity, lower); 279 271 for (NodeIt n(_graph); n != INVALID; ++n) { … … 360 352 /// Run the algorithm. 361 353 /// 354 /// \param partial_augment By default the algorithm performs 355 /// partial augment and relabel operations in the cost scaling 356 /// phases. Set this parameter to \c false for using local push and 357 /// relabel operations instead. 358 /// 362 359 /// \return \c true if a feasible flow can be found. 363 bool run() { 364 return init() && start(); 360 bool run(bool partial_augment = true) { 361 if (partial_augment) { 362 return init() && startPartialAugment(); 363 } else { 364 return init() && startPushRelabel(); 365 } 365 366 } 366 367 … … 434 435 bool init() { 435 436 if (!_valid_supply) return false; 436 437 // Initializing flow and potential maps 437 // The scaling factor 438 _alpha = 8; 439 440 // Initialize flow and potential maps 438 441 if (!_flow) { 439 442 _flow = new FlowMap(_graph); … … 448 451 _res_graph = new ResGraph(_graph, _capacity, *_flow); 449 452 450 // Initializ ingthe scaled cost map and the epsilon parameter453 // Initialize the scaled cost map and the epsilon parameter 451 454 Cost max_cost = 0; 452 455 int node_num = countNodes(_graph); 453 456 for (EdgeIt e(_graph); e != INVALID; ++e) { 454 _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA;457 _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha; 455 458 if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e]; 456 459 } 457 460 _epsilon = max_cost * node_num; 458 461 459 // Find inga feasible flow using Circulation462 // Find a feasible flow using Circulation 460 463 Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap, 461 464 SupplyMap > … … 465 468 } 466 469 467 468 /// Execute the algorithm. 469 bool start() { 470 /// Execute the algorithm performing partial augmentation and 471 /// relabel operations. 472 bool startPartialAugment() { 473 // Paramters for heuristics 474 const int BF_HEURISTIC_EPSILON_BOUND = 1000; 475 const int BF_HEURISTIC_BOUND_FACTOR = 3; 476 // Maximum augment path length 477 const int MAX_PATH_LENGTH = 4; 478 479 // Variables 480 typename Graph::template NodeMap<Edge> pred_edge(_graph); 481 typename Graph::template NodeMap<bool> forward(_graph); 482 typename Graph::template NodeMap<OutEdgeIt> next_out(_graph); 483 typename Graph::template NodeMap<InEdgeIt> next_in(_graph); 484 typename Graph::template NodeMap<bool> next_dir(_graph); 470 485 std::deque<Node> active_nodes; 471 typename Graph::template NodeMap<bool> hyper(_graph, false);486 std::vector<Node> path_nodes; 472 487 473 488 int node_num = countNodes(_graph); 474 for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA&& _epsilon > 1 ?475 1 : _epsilon / ALPHA)489 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 490 1 : _epsilon / _alpha ) 476 491 { 477 // Performing price refinement heuristic using BellmanFord478 // algorithm492 // "Early Termination" heuristic: use BellmanFord algorithm 493 // to check if the current flow is optimal 479 494 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { 480 495 typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap; 481 ShiftCostMap shift_cost(_res_cost, _epsilon);496 ShiftCostMap shift_cost(_res_cost, 1); 482 497 BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost); 483 498 bf.init(0); … … 486 501 for (int i = 0; i < K && !done; ++i) 487 502 done = bf.processNextWeakRound(); 488 if (done) { 489 for (NodeIt n(_graph); n != INVALID; ++n) 490 (*_potential)[n] = bf.dist(n); 491 continue; 492 } 493 } 494 495 // Saturating edges not satisfying the optimality condition 503 if (done) break; 504 } 505 506 // Saturate edges not satisfying the optimality condition 496 507 Capacity delta; 497 508 for (EdgeIt e(_graph); e != INVALID; ++e) { … … 509 520 } 510 521 511 // Find ingactive nodes (i.e. nodes with positive excess)512 for (NodeIt n(_graph); n != INVALID; ++n) 522 // Find active nodes (i.e. nodes with positive excess) 523 for (NodeIt n(_graph); n != INVALID; ++n) { 513 524 if (_excess[n] > 0) active_nodes.push_back(n); 514 515 // Performing push and relabel operations 525 } 526 527 // Initialize the next edge maps 528 for (NodeIt n(_graph); n != INVALID; ++n) { 529 next_out[n] = OutEdgeIt(_graph, n); 530 next_in[n] = InEdgeIt(_graph, n); 531 next_dir[n] = true; 532 } 533 534 // Perform partial augment and relabel operations 516 535 while (active_nodes.size() > 0) { 536 // Select an active node (FIFO selection) 537 if (_excess[active_nodes[0]] <= 0) { 538 active_nodes.pop_front(); 539 continue; 540 } 541 Node start = active_nodes[0]; 542 path_nodes.clear(); 543 path_nodes.push_back(start); 544 545 // Find an augmenting path from the start node 546 Node u, tip = start; 547 LCost min_red_cost; 548 while ( _excess[tip] >= 0 && 549 int(path_nodes.size()) <= MAX_PATH_LENGTH ) 550 { 551 if (next_dir[tip]) { 552 for (OutEdgeIt e = next_out[tip]; e != INVALID; ++e) { 553 if (_capacity[e]  (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { 554 u = _graph.target(e); 555 pred_edge[u] = e; 556 forward[u] = true; 557 next_out[tip] = e; 558 tip = u; 559 path_nodes.push_back(tip); 560 goto next_step; 561 } 562 } 563 next_dir[tip] = false; 564 } 565 for (InEdgeIt e = next_in[tip]; e != INVALID; ++e) { 566 if ((*_flow)[e] > 0 && (*_red_cost)[e] < 0) { 567 u = _graph.source(e); 568 pred_edge[u] = e; 569 forward[u] = false; 570 next_in[tip] = e; 571 tip = u; 572 path_nodes.push_back(tip); 573 goto next_step; 574 } 575 } 576 577 // Relabel tip node 578 min_red_cost = std::numeric_limits<LCost>::max() / 2; 579 for (OutEdgeIt oe(_graph, tip); oe != INVALID; ++oe) { 580 if ( _capacity[oe]  (*_flow)[oe] > 0 && 581 (*_red_cost)[oe] < min_red_cost ) 582 min_red_cost = (*_red_cost)[oe]; 583 } 584 for (InEdgeIt ie(_graph, tip); ie != INVALID; ++ie) { 585 if ((*_flow)[ie] > 0 && (*_red_cost)[ie] < min_red_cost) 586 min_red_cost = (*_red_cost)[ie]; 587 } 588 (*_potential)[tip] = min_red_cost + _epsilon; 589 590 // Reset the next edge maps 591 next_out[tip] = OutEdgeIt(_graph, tip); 592 next_in[tip] = InEdgeIt(_graph, tip); 593 next_dir[tip] = true; 594 595 // Step back 596 if (tip != start) { 597 path_nodes.pop_back(); 598 tip = path_nodes[path_nodes.size()1]; 599 } 600 601 next_step: 602 continue; 603 } 604 605 // Augment along the found path (as much flow as possible) 606 Capacity delta; 607 for (int i = 1; i < int(path_nodes.size()); ++i) { 608 u = path_nodes[i]; 609 delta = forward[u] ? 610 _capacity[pred_edge[u]]  (*_flow)[pred_edge[u]] : 611 (*_flow)[pred_edge[u]]; 612 delta = std::min(delta, _excess[path_nodes[i1]]); 613 (*_flow)[pred_edge[u]] += forward[u] ? delta : delta; 614 _excess[path_nodes[i1]] = delta; 615 _excess[u] += delta; 616 if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u); 617 } 618 } 619 } 620 621 // Compute node potentials for the original costs 622 ResidualCostMap<CostMap> res_cost(_orig_cost); 623 BellmanFord< ResGraph, ResidualCostMap<CostMap> > 624 bf(*_res_graph, res_cost); 625 bf.init(0); bf.start(); 626 for (NodeIt n(_graph); n != INVALID; ++n) 627 (*_potential)[n] = bf.dist(n); 628 629 // Handle nonzero lower bounds 630 if (_lower) { 631 for (EdgeIt e(_graph); e != INVALID; ++e) 632 (*_flow)[e] += (*_lower)[e]; 633 } 634 return true; 635 } 636 637 /// Execute the algorithm performing push and relabel operations. 638 bool startPushRelabel() { 639 // Paramters for heuristics 640 const int BF_HEURISTIC_EPSILON_BOUND = 1000; 641 const int BF_HEURISTIC_BOUND_FACTOR = 3; 642 643 typename Graph::template NodeMap<bool> hyper(_graph, false); 644 typename Graph::template NodeMap<Edge> pred_edge(_graph); 645 typename Graph::template NodeMap<bool> forward(_graph); 646 typename Graph::template NodeMap<OutEdgeIt> next_out(_graph); 647 typename Graph::template NodeMap<InEdgeIt> next_in(_graph); 648 typename Graph::template NodeMap<bool> next_dir(_graph); 649 std::deque<Node> active_nodes; 650 651 int node_num = countNodes(_graph); 652 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 653 1 : _epsilon / _alpha ) 654 { 655 // "Early Termination" heuristic: use BellmanFord algorithm 656 // to check if the current flow is optimal 657 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { 658 typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap; 659 ShiftCostMap shift_cost(_res_cost, 1); 660 BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost); 661 bf.init(0); 662 bool done = false; 663 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); 664 for (int i = 0; i < K && !done; ++i) 665 done = bf.processNextWeakRound(); 666 if (done) break; 667 } 668 669 // Saturate edges not satisfying the optimality condition 670 Capacity delta; 671 for (EdgeIt e(_graph); e != INVALID; ++e) { 672 if (_capacity[e]  (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { 673 delta = _capacity[e]  (*_flow)[e]; 674 _excess[_graph.source(e)] = delta; 675 _excess[_graph.target(e)] += delta; 676 (*_flow)[e] = _capacity[e]; 677 } 678 if ((*_flow)[e] > 0 && (*_red_cost)[e] < 0) { 679 _excess[_graph.target(e)] = (*_flow)[e]; 680 _excess[_graph.source(e)] += (*_flow)[e]; 681 (*_flow)[e] = 0; 682 } 683 } 684 685 // Find active nodes (i.e. nodes with positive excess) 686 for (NodeIt n(_graph); n != INVALID; ++n) { 687 if (_excess[n] > 0) active_nodes.push_back(n); 688 } 689 690 // Initialize the next edge maps 691 for (NodeIt n(_graph); n != INVALID; ++n) { 692 next_out[n] = OutEdgeIt(_graph, n); 693 next_in[n] = InEdgeIt(_graph, n); 694 next_dir[n] = true; 695 } 696 697 // Perform push and relabel operations 698 while (active_nodes.size() > 0) { 699 // Select an active node (FIFO selection) 517 700 Node n = active_nodes[0], t; 518 701 bool relabel_enabled = true; 519 702 520 // Performing push operations if there are admissible edges 521 if (_excess[n] > 0) { 522 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 703 // Perform push operations if there are admissible edges 704 if (_excess[n] > 0 && next_dir[n]) { 705 OutEdgeIt e = next_out[n]; 706 for ( ; e != INVALID; ++e) { 523 707 if (_capacity[e]  (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { 524 delta = _capacity[e]  (*_flow)[e] <= _excess[n] ? 525 _capacity[e]  (*_flow)[e] : _excess[n]; 708 delta = std::min(_capacity[e]  (*_flow)[e], _excess[n]); 526 709 t = _graph.target(e); 527 710 … … 538 721 if (ahead < 0) ahead = 0; 539 722 540 // Push ingflow along the edge723 // Push flow along the edge 541 724 if (ahead < delta) { 542 725 (*_flow)[e] += ahead; … … 558 741 } 559 742 } 560 } 561 562 if (_excess[n] > 0) { 563 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 743 if (e != INVALID) { 744 next_out[n] = e; 745 } else { 746 next_dir[n] = false; 747 } 748 } 749 750 if (_excess[n] > 0 && !next_dir[n]) { 751 InEdgeIt e = next_in[n]; 752 for ( ; e != INVALID; ++e) { 564 753 if ((*_flow)[e] > 0 && (*_red_cost)[e] < 0) { 565 delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];754 delta = std::min((*_flow)[e], _excess[n]); 566 755 t = _graph.source(e); 567 756 … … 578 767 if (ahead < 0) ahead = 0; 579 768 580 // Push ingflow along the edge769 // Push flow along the edge 581 770 if (ahead < delta) { 582 771 (*_flow)[e] = ahead; … … 598 787 } 599 788 } 600 } 601 789 next_in[n] = e; 790 } 791 792 // Relabel the node if it is still active (or hyper) 602 793 if (relabel_enabled && (_excess[n] > 0  hyper[n])) { 603 // Performing relabel operation if the node is still active 604 LCost min_red_cost = std::numeric_limits<LCost>::max(); 794 LCost min_red_cost = std::numeric_limits<LCost>::max() / 2; 605 795 for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) { 606 796 if ( _capacity[oe]  (*_flow)[oe] > 0 && … … 614 804 (*_potential)[n] = min_red_cost + _epsilon; 615 805 hyper[n] = false; 616 } 617 618 // Removing active nodes with nonpositive excess 806 807 // Reset the next edge maps 808 next_out[n] = OutEdgeIt(_graph, n); 809 next_in[n] = InEdgeIt(_graph, n); 810 next_dir[n] = true; 811 } 812 813 // Remove nodes that are not active nor hyper 619 814 while ( active_nodes.size() > 0 && 620 815 _excess[active_nodes[0]] <= 0 && … … 625 820 } 626 821 627 // Comput ingnode potentials for the original costs822 // Compute node potentials for the original costs 628 823 ResidualCostMap<CostMap> res_cost(_orig_cost); 629 824 BellmanFord< ResGraph, ResidualCostMap<CostMap> > … … 633 828 (*_potential)[n] = bf.dist(n); 634 829 635 // Handl ingnonzero lower bounds830 // Handle nonzero lower bounds 636 831 if (_lower) { 637 832 for (EdgeIt e(_graph); e != INVALID; ++e)
Note: See TracChangeset
for help on using the changeset viewer.