Changeset 2634:e98bbe64cca4 in lemon0.x
 Timestamp:
 02/06/09 22:52:34 (16 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@3519
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/network_simplex.h
r2630 r2634 29 29 #include <algorithm> 30 30 31 #include <lemon/graph_adaptor.h>32 31 #include <lemon/graph_utils.h> 33 #include <lemon/smart_graph.h>34 32 #include <lemon/math.h> 35 33 … … 73 71 class NetworkSimplex 74 72 { 73 GRAPH_TYPEDEFS(typename Graph); 74 75 75 typedef typename CapacityMap::Value Capacity; 76 76 typedef typename CostMap::Value Cost; 77 77 typedef typename SupplyMap::Value Supply; 78 78 79 typedef SmartGraph SGraph;80 GRAPH_TYPEDEFS(typename SGraph);81 82 typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;83 typedef typename SGraph::template EdgeMap<Cost> SCostMap;84 typedef typename SGraph::template NodeMap<Supply> SSupplyMap;85 typedef typename SGraph::template NodeMap<Cost> SPotentialMap;86 87 typedef typename SGraph::template NodeMap<int> IntNodeMap;88 typedef typename SGraph::template NodeMap<bool> BoolNodeMap;89 typedef typename SGraph::template NodeMap<Node> NodeNodeMap;90 typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;91 typedef typename SGraph::template EdgeMap<int> IntEdgeMap;92 typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;93 94 typedef typename Graph::template NodeMap<Node> NodeRefMap;95 typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;96 97 79 typedef std::vector<Edge> EdgeVector; 80 typedef std::vector<Node> NodeVector; 81 typedef std::vector<int> IntVector; 82 typedef std::vector<bool> BoolVector; 83 typedef std::vector<Capacity> CapacityVector; 84 typedef std::vector<Cost> CostVector; 85 typedef std::vector<Supply> SupplyVector; 86 87 typedef typename Graph::template NodeMap<int> IntNodeMap; 98 88 99 89 public: … … 117 107 private: 118 108 119 /// \brief Map adaptor class for handling reduced edge costs.120 ///121 /// Map adaptor class for handling reduced edge costs.122 class ReducedCostMap : public MapBase<Edge, Cost>123 {124 private:125 126 const SGraph &_gr;127 const SCostMap &_cost_map;128 const SPotentialMap &_pot_map;129 130 public:131 132 ///\e133 ReducedCostMap( const SGraph &gr,134 const SCostMap &cost_map,135 const SPotentialMap &pot_map ) :136 _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}137 138 ///\e139 Cost operator[](const Edge &e) const {140 return _cost_map[e] + _pot_map[_gr.source(e)]141  _pot_map[_gr.target(e)];142 }143 144 }; //class ReducedCostMap145 146 private:147 148 109 /// \brief Implementation of the "First Eligible" pivot rule for the 149 110 /// \ref NetworkSimplex "network simplex" algorithm. … … 158 119 159 120 // References to the NetworkSimplex class 160 NetworkSimplex &_ns; 161 EdgeVector &_edges; 162 121 const EdgeVector &_edge; 122 const IntVector &_source; 123 const IntVector &_target; 124 const CostVector &_cost; 125 const IntVector &_state; 126 const CostVector &_pi; 127 int &_in_edge; 128 int _edge_num; 129 130 // Pivot rule data 163 131 int _next_edge; 164 132 … … 166 134 167 135 /// Constructor 168 FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) : 169 _ns(ns), _edges(edges), _next_edge(0) {} 136 FirstEligiblePivotRule(NetworkSimplex &ns) : 137 _edge(ns._edge), _source(ns._source), _target(ns._target), 138 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 139 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0) 140 {} 170 141 171 142 /// Find next entering edge 172 143 bool findEnteringEdge() { 173 Edge e;174 for (int i = _next_edge; i < int(_edges.size()); ++i) {175 e = _edges[i];176 if ( _ns._state[e] * _ns._red_cost[e]< 0) {177 _ ns._in_edge = e;178 _next_edge = i+ 1;144 Cost c; 145 for (int e = _next_edge; e < _edge_num; ++e) { 146 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 147 if (c < 0) { 148 _in_edge = e; 149 _next_edge = e + 1; 179 150 return true; 180 151 } 181 152 } 182 for (int i = 0; i < _next_edge; ++i) {183 e = _edges[i];184 if ( _ns._state[e] * _ns._red_cost[e]< 0) {185 _ ns._in_edge = e;186 _next_edge = i+ 1;153 for (int e = 0; e < _next_edge; ++e) { 154 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 155 if (c < 0) { 156 _in_edge = e; 157 _next_edge = e + 1; 187 158 return true; 188 159 } … … 190 161 return false; 191 162 } 163 192 164 }; //class FirstEligiblePivotRule 165 193 166 194 167 /// \brief Implementation of the "Best Eligible" pivot rule for the … … 204 177 205 178 // References to the NetworkSimplex class 206 NetworkSimplex &_ns; 207 EdgeVector &_edges; 179 const EdgeVector &_edge; 180 const IntVector &_source; 181 const IntVector &_target; 182 const CostVector &_cost; 183 const IntVector &_state; 184 const CostVector &_pi; 185 int &_in_edge; 186 int _edge_num; 208 187 209 188 public: 210 189 211 190 /// Constructor 212 BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) : 213 _ns(ns), _edges(edges) {} 191 BestEligiblePivotRule(NetworkSimplex &ns) : 192 _edge(ns._edge), _source(ns._source), _target(ns._target), 193 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 194 _in_edge(ns._in_edge), _edge_num(ns._edge_num) 195 {} 214 196 215 197 /// Find next entering edge 216 198 bool findEnteringEdge() { 217 Cost min = 0; 218 Edge e; 219 for (int i = 0; i < int(_edges.size()); ++i) { 220 e = _edges[i]; 221 if (_ns._state[e] * _ns._red_cost[e] < min) { 222 min = _ns._state[e] * _ns._red_cost[e]; 223 _ns._in_edge = e; 199 Cost c, min = 0; 200 for (int e = 0; e < _edge_num; ++e) { 201 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 202 if (c < min) { 203 min = c; 204 _in_edge = e; 224 205 } 225 206 } 226 207 return min < 0; 227 208 } 209 228 210 }; //class BestEligiblePivotRule 211 229 212 230 213 /// \brief Implementation of the "Block Search" pivot rule for the … … 240 223 241 224 // References to the NetworkSimplex class 242 NetworkSimplex &_ns; 243 EdgeVector &_edges; 244 225 const EdgeVector &_edge; 226 const IntVector &_source; 227 const IntVector &_target; 228 const CostVector &_cost; 229 const IntVector &_state; 230 const CostVector &_pi; 231 int &_in_edge; 232 int _edge_num; 233 234 // Pivot rule data 245 235 int _block_size; 246 int _next_edge , _min_edge;236 int _next_edge; 247 237 248 238 public: 249 239 250 240 /// Constructor 251 BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) : 252 _ns(ns), _edges(edges), _next_edge(0), _min_edge(0) 241 BlockSearchPivotRule(NetworkSimplex &ns) : 242 _edge(ns._edge), _source(ns._source), _target(ns._target), 243 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 244 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0) 253 245 { 254 246 // The main parameters of the pivot rule … … 256 248 const int MIN_BLOCK_SIZE = 10; 257 249 258 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge s.size())),250 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)), 259 251 MIN_BLOCK_SIZE ); 260 252 } … … 262 254 /// Find next entering edge 263 255 bool findEnteringEdge() { 264 Cost curr, min = 0; 265 Edge e; 256 Cost c, min = 0; 266 257 int cnt = _block_size; 267 int i;268 for ( i = _next_edge; i < int(_edges.size()); ++i) {269 e = _edges[i];270 if ( (curr = _ns._state[e] * _ns._red_cost[e])< min) {271 min = c urr;272 _min_edge = i;258 int e, min_edge = _next_edge; 259 for (e = _next_edge; e < _edge_num; ++e) { 260 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 261 if (c < min) { 262 min = c; 263 min_edge = e; 273 264 } 274 265 if (cnt == 0) { … … 278 269 } 279 270 if (min == 0  cnt > 0) { 280 for ( i = 0; i < _next_edge; ++i) {281 e = _edges[i];282 if ( (curr = _ns._state[e] * _ns._red_cost[e])< min) {283 min = c urr;284 _min_edge = i;271 for (e = 0; e < _next_edge; ++e) { 272 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 273 if (c < min) { 274 min = c; 275 min_edge = e; 285 276 } 286 277 if (cnt == 0) { … … 291 282 } 292 283 if (min >= 0) return false; 293 _ ns._in_edge = _edges[_min_edge];294 _next_edge = i;284 _in_edge = min_edge; 285 _next_edge = e; 295 286 return true; 296 287 } 288 297 289 }; //class BlockSearchPivotRule 290 298 291 299 292 /// \brief Implementation of the "Candidate List" pivot rule for the … … 309 302 310 303 // References to the NetworkSimplex class 311 NetworkSimplex &_ns; 312 EdgeVector &_edges; 313 314 EdgeVector _candidates; 304 const EdgeVector &_edge; 305 const IntVector &_source; 306 const IntVector &_target; 307 const CostVector &_cost; 308 const IntVector &_state; 309 const CostVector &_pi; 310 int &_in_edge; 311 int _edge_num; 312 313 // Pivot rule data 314 IntVector _candidates; 315 315 int _list_length, _minor_limit; 316 316 int _curr_length, _minor_count; 317 int _next_edge , _min_edge;317 int _next_edge; 318 318 319 319 public: 320 320 321 321 /// Constructor 322 CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) : 323 _ns(ns), _edges(edges), _next_edge(0), _min_edge(0) 322 CandidateListPivotRule(NetworkSimplex &ns) : 323 _edge(ns._edge), _source(ns._source), _target(ns._target), 324 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 325 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0) 324 326 { 325 327 // The main parameters of the pivot rule … … 329 331 const int MIN_MINOR_LIMIT = 3; 330 332 331 _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge s.size())),333 _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)), 332 334 MIN_LIST_LENGTH ); 333 335 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), … … 339 341 /// Find next entering edge 340 342 bool findEnteringEdge() { 341 Cost min, curr; 343 Cost min, c; 344 int e, min_edge = _next_edge; 342 345 if (_curr_length > 0 && _minor_count < _minor_limit) { 343 346 // Minor iteration: select the best eligible edge from the 344 347 // current candidate list 345 348 ++_minor_count; 346 Edge e;347 349 min = 0; 348 350 for (int i = 0; i < _curr_length; ++i) { 349 351 e = _candidates[i]; 350 c urr = _ns._state[e] * _ns._red_cost[e];351 if (c urr< min) {352 min = c urr;353 _ns._in_edge = e;352 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 353 if (c < min) { 354 min = c; 355 min_edge = e; 354 356 } 355 if (c urr>= 0) {357 if (c >= 0) { 356 358 _candidates[i] = _candidates[_curr_length]; 357 359 } 358 360 } 359 if (min < 0) return true; 361 if (min < 0) { 362 _in_edge = min_edge; 363 return true; 364 } 360 365 } 361 366 362 367 // Major iteration: build a new candidate list 363 Edge e;364 368 min = 0; 365 369 _curr_length = 0; 366 int i; 367 for (i = _next_edge; i < int(_edges.size()); ++i) { 368 e = _edges[i]; 369 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { 370 for (e = _next_edge; e < _edge_num; ++e) { 371 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 372 if (c < 0) { 370 373 _candidates[_curr_length++] = e; 371 if (c urr< min) {372 min = c urr;373 _min_edge = i;374 if (c < min) { 375 min = c; 376 min_edge = e; 374 377 } 375 378 if (_curr_length == _list_length) break; … … 377 380 } 378 381 if (_curr_length < _list_length) { 379 for ( i = 0; i < _next_edge; ++i) {380 e = _edges[i];381 if ( (curr = _ns._state[e] * _ns._red_cost[e])< 0) {382 for (e = 0; e < _next_edge; ++e) { 383 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 384 if (c < 0) { 382 385 _candidates[_curr_length++] = e; 383 if (c urr< min) {384 min = c urr;385 _min_edge = i;386 if (c < min) { 387 min = c; 388 min_edge = e; 386 389 } 387 390 if (_curr_length == _list_length) break; … … 391 394 if (_curr_length == 0) return false; 392 395 _minor_count = 1; 393 _ ns._in_edge = _edges[_min_edge];394 _next_edge = i;396 _in_edge = min_edge; 397 _next_edge = e; 395 398 return true; 396 399 } 400 397 401 }; //class CandidateListPivotRule 402 398 403 399 404 /// \brief Implementation of the "Altering Candidate List" pivot rule … … 409 414 410 415 // References to the NetworkSimplex class 411 NetworkSimplex &_ns; 412 EdgeVector &_edges; 413 414 EdgeVector _candidates; 415 SCostMap _cand_cost; 416 const EdgeVector &_edge; 417 const IntVector &_source; 418 const IntVector &_target; 419 const CostVector &_cost; 420 const IntVector &_state; 421 const CostVector &_pi; 422 int &_in_edge; 423 int _edge_num; 424 416 425 int _block_size, _head_length, _curr_length; 417 426 int _next_edge; 427 IntVector _candidates; 428 CostVector _cand_cost; 418 429 419 430 // Functor class to compare edges during sort of the candidate list … … 421 432 { 422 433 private: 423 const SCostMap&_map;434 const CostVector &_map; 424 435 public: 425 SortFunc(const SCostMap&map) : _map(map) {}426 bool operator()( const Edge &e1, const Edge &e2) {427 return _map[ e1] > _map[e2];436 SortFunc(const CostVector &map) : _map(map) {} 437 bool operator()(int left, int right) { 438 return _map[left] > _map[right]; 428 439 } 429 440 }; … … 434 445 435 446 /// Constructor 436 AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) : 437 _ns(ns), _edges(edges), _cand_cost(_ns._graph), 438 _next_edge(0), _sort_func(_cand_cost) 447 AlteringListPivotRule(NetworkSimplex &ns) : 448 _edge(ns._edge), _source(ns._source), _target(ns._target), 449 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 450 _in_edge(ns._in_edge), _edge_num(ns._edge_num), 451 _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost) 439 452 { 440 453 // The main parameters of the pivot rule … … 444 457 const int MIN_HEAD_LENGTH = 3; 445 458 446 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge s.size())),459 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)), 447 460 MIN_BLOCK_SIZE ); 448 461 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), … … 455 468 bool findEnteringEdge() { 456 469 // Check the current candidate list 457 Edge e; 458 for (int ix = 0; ix < _curr_length; ++ix) { 459 e = _candidates[ix]; 460 if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) { 461 _candidates[ix] = _candidates[_curr_length]; 470 int e; 471 for (int i = 0; i < _curr_length; ++i) { 472 e = _candidates[i]; 473 _cand_cost[e] = _state[e] * 474 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 475 if (_cand_cost[e] >= 0) { 476 _candidates[i] = _candidates[_curr_length]; 462 477 } 463 478 } … … 467 482 int last_edge = 0; 468 483 int limit = _head_length; 469 for (int i = _next_edge; i < int(_edges.size()); ++i) { 470 e = _edges[i]; 471 if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) { 484 485 for (int e = _next_edge; e < _edge_num; ++e) { 486 _cand_cost[e] = _state[e] * 487 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 488 if (_cand_cost[e] < 0) { 472 489 _candidates[_curr_length++] = e; 473 last_edge = i;490 last_edge = e; 474 491 } 475 492 if (cnt == 0) { … … 480 497 } 481 498 if (_curr_length <= limit) { 482 for (int i = 0; i < _next_edge; ++i) { 483 e = _edges[i]; 484 if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) { 499 for (int e = 0; e < _next_edge; ++e) { 500 _cand_cost[e] = _state[e] * 501 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 502 if (_cand_cost[e] < 0) { 485 503 _candidates[_curr_length++] = e; 486 last_edge = i;504 last_edge = e; 487 505 } 488 506 if (cnt == 0) { … … 501 519 502 520 // Pop the first element of the heap 503 _ ns._in_edge = _candidates[0];521 _in_edge = _candidates[0]; 504 522 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, 505 523 _sort_func ); … … 507 525 return true; 508 526 } 527 509 528 }; //class AlteringListPivotRule 510 529 … … 520 539 private: 521 540 522 // The directed graph the algorithm runs on523 SGraph _graph;524 541 // The original graph 525 const Graph &_ graph_ref;542 const Graph &_orig_graph; 526 543 // The original lower bound map 527 const LowerMap *_lower; 528 // The capacity map 529 SCapacityMap _capacity; 530 // The cost map 531 SCostMap _cost; 532 // The supply map 533 SSupplyMap _supply; 534 bool _valid_supply; 535 536 // Edge map of the current flow 537 SCapacityMap _flow; 538 // Node map of the current potentials 539 SPotentialMap _potential; 540 541 // The depth node map of the spanning tree structure 542 IntNodeMap _depth; 543 // The parent node map of the spanning tree structure 544 NodeNodeMap _parent; 545 // The pred_edge node map of the spanning tree structure 546 EdgeNodeMap _pred_edge; 547 // The thread node map of the spanning tree structure 548 NodeNodeMap _thread; 549 // The forward node map of the spanning tree structure 550 BoolNodeMap _forward; 551 // The state edge map 552 IntEdgeMap _state; 553 // The artificial root node of the spanning tree structure 554 Node _root; 555 556 // The reduced cost map 557 ReducedCostMap _red_cost; 558 559 // The nonartifical edges 560 EdgeVector _edges; 561 562 // Members for handling the original graph 544 const LowerMap *_orig_lower; 545 // The original capacity map 546 const CapacityMap &_orig_cap; 547 // The original cost map 548 const CostMap &_orig_cost; 549 // The original supply map 550 const SupplyMap *_orig_supply; 551 // The source node (if no supply map was given) 552 Node _orig_source; 553 // The target node (if no supply map was given) 554 Node _orig_target; 555 // The flow value (if no supply map was given) 556 Capacity _orig_flow_value; 557 558 // The flow result map 563 559 FlowMap *_flow_result; 560 // The potential result map 564 561 PotentialMap *_potential_result; 562 // Indicate if the flow result map is local 565 563 bool _local_flow; 564 // Indicate if the potential result map is local 566 565 bool _local_potential; 567 NodeRefMap _node_ref; 568 EdgeRefMap _edge_ref; 566 567 // The edge references 568 EdgeVector _edge; 569 // The node references 570 NodeVector _node; 571 // The node ids 572 IntNodeMap _node_id; 573 // The source nodes of the edges 574 IntVector _source; 575 // The target nodess of the edges 576 IntVector _target; 577 578 // The (modified) capacity vector 579 CapacityVector _cap; 580 // The cost vector 581 CostVector _cost; 582 // The (modified) supply vector 583 CostVector _supply; 584 // The current flow vector 585 CapacityVector _flow; 586 // The current potential vector 587 CostVector _pi; 588 589 // The number of nodes in the original graph 590 int _node_num; 591 // The number of edges in the original graph 592 int _edge_num; 593 594 // The depth vector of the spanning tree structure 595 IntVector _depth; 596 // The parent vector of the spanning tree structure 597 IntVector _parent; 598 // The pred_edge vector of the spanning tree structure 599 IntVector _pred; 600 // The thread vector of the spanning tree structure 601 IntVector _thread; 602 // The forward vector of the spanning tree structure 603 BoolVector _forward; 604 // The state vector 605 IntVector _state; 606 // The root node 607 int _root; 569 608 570 609 // The entering edge of the current pivot iteration 571 Edge_in_edge;610 int _in_edge; 572 611 573 612 // Temporary nodes used in the current pivot iteration 574 Node join, u_in, v_in, u_out, v_out; 575 Node right, first, second, last; 576 Node stem, par_stem, new_stem; 613 int join, u_in, v_in, u_out, v_out; 614 int right, first, second, last; 615 int stem, par_stem, new_stem; 616 577 617 // The maximum augment amount along the found cycle in the current 578 618 // pivot iteration 579 619 Capacity delta; 580 620 581 public 621 public: 582 622 583 623 /// \brief General constructor (with lower bounds). … … 595 635 const CostMap &cost, 596 636 const SupplyMap &supply ) : 597 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), 598 _cost(_graph), _supply(_graph), _flow(_graph), 599 _potential(_graph), _depth(_graph), _parent(_graph), 600 _pred_edge(_graph), _thread(_graph), _forward(_graph), 601 _state(_graph), _red_cost(_graph, _cost, _potential), 637 _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity), 638 _orig_cost(cost), _orig_supply(&supply), 602 639 _flow_result(NULL), _potential_result(NULL), 603 640 _local_flow(false), _local_potential(false), 604 _node_ref(graph), _edge_ref(graph) 605 { 606 // Check the sum of supply values 607 Supply sum = 0; 608 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) 609 sum += supply[n]; 610 if (!(_valid_supply = sum == 0)) return; 611 612 // Copy _graph_ref to _graph 613 _graph.reserveNode(countNodes(_graph_ref) + 1); 614 _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); 615 copyGraph(_graph, _graph_ref) 616 .edgeMap(_capacity, capacity) 617 .edgeMap(_cost, cost) 618 .nodeMap(_supply, supply) 619 .nodeRef(_node_ref) 620 .edgeRef(_edge_ref) 621 .run(); 622 623 // Remove nonzero lower bounds 624 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { 625 if (lower[e] != 0) { 626 _capacity[_edge_ref[e]] = capacity[e]  lower[e]; 627 _supply[_node_ref[_graph_ref.source(e)]] = lower[e]; 628 _supply[_node_ref[_graph_ref.target(e)]] += lower[e]; 629 } 630 } 631 } 641 _node_id(graph) 642 {} 632 643 633 644 /// \brief General constructor (without lower bounds). … … 643 654 const CostMap &cost, 644 655 const SupplyMap &supply ) : 645 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), 646 _cost(_graph), _supply(_graph), _flow(_graph), 647 _potential(_graph), _depth(_graph), _parent(_graph), 648 _pred_edge(_graph), _thread(_graph), _forward(_graph), 649 _state(_graph), _red_cost(_graph, _cost, _potential), 656 _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity), 657 _orig_cost(cost), _orig_supply(&supply), 650 658 _flow_result(NULL), _potential_result(NULL), 651 659 _local_flow(false), _local_potential(false), 652 _node_ref(graph), _edge_ref(graph) 653 { 654 // Check the sum of supply values 655 Supply sum = 0; 656 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) 657 sum += supply[n]; 658 if (!(_valid_supply = sum == 0)) return; 659 660 // Copy _graph_ref to _graph 661 _graph.reserveNode(countNodes(_graph_ref) + 1); 662 _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); 663 copyGraph(_graph, _graph_ref) 664 .edgeMap(_capacity, capacity) 665 .edgeMap(_cost, cost) 666 .nodeMap(_supply, supply) 667 .nodeRef(_node_ref) 668 .edgeRef(_edge_ref) 669 .run(); 670 } 660 _node_id(graph) 661 {} 671 662 672 663 /// \brief Simple constructor (with lower bounds). … … 686 677 const CapacityMap &capacity, 687 678 const CostMap &cost, 688 typename Graph::Node s, 689 typename Graph::Node t, 690 typename SupplyMap::Value flow_value ) : 691 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), 692 _cost(_graph), _supply(_graph, 0), _flow(_graph), 693 _potential(_graph), _depth(_graph), _parent(_graph), 694 _pred_edge(_graph), _thread(_graph), _forward(_graph), 695 _state(_graph), _red_cost(_graph, _cost, _potential), 679 Node s, Node t, 680 Capacity flow_value ) : 681 _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity), 682 _orig_cost(cost), _orig_supply(NULL), 683 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), 696 684 _flow_result(NULL), _potential_result(NULL), 697 685 _local_flow(false), _local_potential(false), 698 _node_ref(graph), _edge_ref(graph) 699 { 700 // Copy _graph_ref to graph 701 _graph.reserveNode(countNodes(_graph_ref) + 1); 702 _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); 703 copyGraph(_graph, _graph_ref) 704 .edgeMap(_capacity, capacity) 705 .edgeMap(_cost, cost) 706 .nodeRef(_node_ref) 707 .edgeRef(_edge_ref) 708 .run(); 709 710 // Remove nonzero lower bounds 711 _supply[_node_ref[s]] = flow_value; 712 _supply[_node_ref[t]] = flow_value; 713 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { 714 if (lower[e] != 0) { 715 _capacity[_edge_ref[e]] = capacity[e]  lower[e]; 716 _supply[_node_ref[_graph_ref.source(e)]] = lower[e]; 717 _supply[_node_ref[_graph_ref.target(e)]] += lower[e]; 718 } 719 } 720 _valid_supply = true; 721 } 686 _node_id(graph) 687 {} 722 688 723 689 /// \brief Simple constructor (without lower bounds). … … 735 701 const CapacityMap &capacity, 736 702 const CostMap &cost, 737 typename Graph::Node s, 738 typename Graph::Node t, 739 typename SupplyMap::Value flow_value ) : 740 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), 741 _cost(_graph), _supply(_graph, 0), _flow(_graph), 742 _potential(_graph), _depth(_graph), _parent(_graph), 743 _pred_edge(_graph), _thread(_graph), _forward(_graph), 744 _state(_graph), _red_cost(_graph, _cost, _potential), 703 Node s, Node t, 704 Capacity flow_value ) : 705 _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity), 706 _orig_cost(cost), _orig_supply(NULL), 707 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), 745 708 _flow_result(NULL), _potential_result(NULL), 746 709 _local_flow(false), _local_potential(false), 747 _node_ref(graph), _edge_ref(graph) 748 { 749 // Copy _graph_ref to graph 750 _graph.reserveNode(countNodes(_graph_ref) + 1); 751 _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); 752 copyGraph(_graph, _graph_ref) 753 .edgeMap(_capacity, capacity) 754 .edgeMap(_cost, cost) 755 .nodeRef(_node_ref) 756 .edgeRef(_edge_ref) 757 .run(); 758 _supply[_node_ref[s]] = flow_value; 759 _supply[_node_ref[t]] = flow_value; 760 _valid_supply = true; 761 } 710 _node_id(graph) 711 {} 762 712 763 713 /// Destructor. … … 895 845 Cost totalCost() const { 896 846 Cost c = 0; 897 for ( typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)898 c += (*_flow_result)[e] * _ cost[_edge_ref[e]];847 for (EdgeIt e(_orig_graph); e != INVALID; ++e) 848 c += (*_flow_result)[e] * _orig_cost[e]; 899 849 return c; 900 850 } … … 904 854 private: 905 855 906 // Extend the underlying graph and initialize all the node and edge maps856 // Initialize internal data structures 907 857 bool init() { 908 if (!_valid_supply) return false;909 910 858 // Initialize result maps 911 859 if (!_flow_result) { 912 _flow_result = new FlowMap(_ graph_ref);860 _flow_result = new FlowMap(_orig_graph); 913 861 _local_flow = true; 914 862 } 915 863 if (!_potential_result) { 916 _potential_result = new PotentialMap(_ graph_ref);864 _potential_result = new PotentialMap(_orig_graph); 917 865 _local_potential = true; 918 866 } 919 920 // Initialize the edge vector and the edge maps 921 _edges.reserve(countEdges(_graph)); 922 for (EdgeIt e(_graph); e != INVALID; ++e) { 923 _edges.push_back(e); 924 _flow[e] = 0; 925 _state[e] = STATE_LOWER; 926 } 927 928 // Add an artificial root node to the graph 929 _root = _graph.addNode(); 930 _parent[_root] = INVALID; 931 _pred_edge[_root] = INVALID; 867 868 // Initialize vectors 869 _node_num = countNodes(_orig_graph); 870 _edge_num = countEdges(_orig_graph); 871 int all_node_num = _node_num + 1; 872 int all_edge_num = _edge_num + _node_num; 873 874 _edge.resize(_edge_num); 875 _node.reserve(_node_num); 876 _source.resize(all_edge_num); 877 _target.resize(all_edge_num); 878 879 _cap.resize(all_edge_num); 880 _cost.resize(all_edge_num); 881 _supply.resize(all_node_num); 882 _flow.resize(all_edge_num, 0); 883 _pi.resize(all_node_num, 0); 884 885 _depth.resize(all_node_num); 886 _parent.resize(all_node_num); 887 _pred.resize(all_node_num); 888 _thread.resize(all_node_num); 889 _forward.resize(all_node_num); 890 _state.resize(all_edge_num, STATE_LOWER); 891 892 // Initialize node related data 893 bool valid_supply = true; 894 if (_orig_supply) { 895 Supply sum = 0; 896 int i = 0; 897 for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) { 898 _node.push_back(n); 899 _node_id[n] = i; 900 _supply[i] = (*_orig_supply)[n]; 901 sum += _supply[i]; 902 } 903 valid_supply = (sum == 0); 904 } else { 905 int i = 0; 906 for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) { 907 _node.push_back(n); 908 _node_id[n] = i; 909 _supply[i] = 0; 910 } 911 _supply[_node_id[_orig_source]] = _orig_flow_value; 912 _supply[_node_id[_orig_target]] = _orig_flow_value; 913 } 914 if (!valid_supply) return false; 915 916 // Set data for the artificial root node 917 _root = _node_num; 932 918 _depth[_root] = 0; 919 _parent[_root] = 1; 920 _pred[_root] = 1; 921 _thread[_root] = 0; 933 922 _supply[_root] = 0; 934 _potential[_root] = 0; 935 936 // Add artificial edges to the graph and initialize the node maps 937 // of the spanning tree data structure 938 Node last = _root; 939 Edge e; 923 _pi[_root] = 0; 924 925 // Store the edges in a mixed order 926 int k = std::max(int(sqrt(_edge_num)), 10); 927 int i = 0; 928 for (EdgeIt e(_orig_graph); e != INVALID; ++e) { 929 _edge[i] = e; 930 if ((i += k) >= _edge_num) i = (i % k) + 1; 931 } 932 933 // Initialize edge maps 934 for (int i = 0; i != _edge_num; ++i) { 935 Edge e = _edge[i]; 936 _source[i] = _node_id[_orig_graph.source(e)]; 937 _target[i] = _node_id[_orig_graph.target(e)]; 938 _cost[i] = _orig_cost[e]; 939 _cap[i] = _orig_cap[e]; 940 } 941 942 // Remove nonzero lower bounds 943 if (_orig_lower) { 944 for (int i = 0; i != _edge_num; ++i) { 945 Capacity c = (*_orig_lower)[_edge[i]]; 946 if (c != 0) { 947 _cap[i] = c; 948 _supply[_source[i]] = c; 949 _supply[_target[i]] += c; 950 } 951 } 952 } 953 954 // Add artificial edges and initialize the spanning tree data structure 940 955 Cost max_cost = std::numeric_limits<Cost>::max() / 4; 941 for (NodeIt u(_graph); u != INVALID; ++u) {942 if (u == _root) continue;943 _thread[ last] = u;944 last = u;956 Capacity max_cap = std::numeric_limits<Capacity>::max(); 957 for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) { 958 _thread[u] = u + 1; 959 _depth[u] = 1; 945 960 _parent[u] = _root; 946 _ depth[u] = 1;961 _pred[u] = e; 947 962 if (_supply[u] >= 0) { 948 e = _graph.addEdge(u, _root);949 963 _flow[e] = _supply[u]; 950 964 _forward[u] = true; 951 _p otential[u] = max_cost;965 _pi[u] = max_cost; 952 966 } else { 953 e = _graph.addEdge(_root, u);954 967 _flow[e] = _supply[u]; 955 968 _forward[u] = false; 956 _p otential[u] = max_cost;969 _pi[u] = max_cost; 957 970 } 958 971 _cost[e] = max_cost; 959 _cap acity[e] = std::numeric_limits<Capacity>::max();972 _cap[e] = max_cap; 960 973 _state[e] = STATE_TREE; 961 _pred_edge[u] = e; 962 } 963 _thread[last] = _root; 974 } 964 975 965 976 return true; … … 968 979 // Find the join node 969 980 void findJoinNode() { 970 Node u = _graph.source(_in_edge); 971 Node v = _graph.target(_in_edge); 981 int u = _source[_in_edge]; 982 int v = _target[_in_edge]; 983 while (_depth[u] > _depth[v]) u = _parent[u]; 984 while (_depth[v] > _depth[u]) v = _parent[v]; 972 985 while (u != v) { 973 if (_depth[u] == _depth[v]) { 974 u = _parent[u]; 975 v = _parent[v]; 976 } 977 else if (_depth[u] > _depth[v]) u = _parent[u]; 978 else v = _parent[v]; 986 u = _parent[u]; 987 v = _parent[v]; 979 988 } 980 989 join = u; 981 990 } 982 991 983 // Find the leaving edge of the cycle and returns true if the 992 // Find the leaving edge of the cycle and returns true if the 984 993 // leaving edge is not the same as the entering edge 985 994 bool findLeavingEdge() { … … 987 996 // of the cycle 988 997 if (_state[_in_edge] == STATE_LOWER) { 989 first = _ graph.source(_in_edge);990 second = _ graph.target(_in_edge);998 first = _source[_in_edge]; 999 second = _target[_in_edge]; 991 1000 } else { 992 first = _ graph.target(_in_edge);993 second = _ graph.source(_in_edge);994 } 995 delta = _cap acity[_in_edge];996 bool result = false;1001 first = _target[_in_edge]; 1002 second = _source[_in_edge]; 1003 } 1004 delta = _cap[_in_edge]; 1005 int result = 0; 997 1006 Capacity d; 998 Edgee;1007 int e; 999 1008 1000 1009 // Search the cycle along the path form the first node to the root 1001 for ( Nodeu = first; u != join; u = _parent[u]) {1002 e = _pred _edge[u];1003 d = _forward[u] ? _flow[e] : _cap acity[e]  _flow[e];1010 for (int u = first; u != join; u = _parent[u]) { 1011 e = _pred[u]; 1012 d = _forward[u] ? _flow[e] : _cap[e]  _flow[e]; 1004 1013 if (d < delta) { 1005 1014 delta = d; 1006 1015 u_out = u; 1007 u_in = first; 1008 v_in = second; 1009 result = true; 1016 result = 1; 1010 1017 } 1011 1018 } 1012 1019 // Search the cycle along the path form the second node to the root 1013 for ( Nodeu = second; u != join; u = _parent[u]) {1014 e = _pred _edge[u];1015 d = _forward[u] ? _cap acity[e]  _flow[e] : _flow[e];1020 for (int u = second; u != join; u = _parent[u]) { 1021 e = _pred[u]; 1022 d = _forward[u] ? _cap[e]  _flow[e] : _flow[e]; 1016 1023 if (d <= delta) { 1017 1024 delta = d; 1018 1025 u_out = u; 1019 u_in = second; 1020 v_in = first; 1021 result = true; 1022 } 1023 } 1024 return result; 1025 } 1026 1027 // Change _flow and _state edge maps 1028 void changeFlows(bool change) { 1026 result = 2; 1027 } 1028 } 1029 1030 if (result == 1) { 1031 u_in = first; 1032 v_in = second; 1033 } else { 1034 u_in = second; 1035 v_in = first; 1036 } 1037 return result != 0; 1038 } 1039 1040 // Change _flow and _state vectors 1041 void changeFlow(bool change) { 1029 1042 // Augment along the cycle 1030 1043 if (delta > 0) { 1031 1044 Capacity val = _state[_in_edge] * delta; 1032 1045 _flow[_in_edge] += val; 1033 for ( Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {1034 _flow[_pred _edge[u]] += _forward[u] ? val : val;1035 } 1036 for ( Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {1037 _flow[_pred _edge[u]] += _forward[u] ? val : val;1046 for (int u = _source[_in_edge]; u != join; u = _parent[u]) { 1047 _flow[_pred[u]] += _forward[u] ? val : val; 1048 } 1049 for (int u = _target[_in_edge]; u != join; u = _parent[u]) { 1050 _flow[_pred[u]] += _forward[u] ? val : val; 1038 1051 } 1039 1052 } … … 1041 1054 if (change) { 1042 1055 _state[_in_edge] = STATE_TREE; 1043 _state[_pred _edge[u_out]] =1044 (_flow[_pred _edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;1056 _state[_pred[u_out]] = 1057 (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; 1045 1058 } else { 1046 1059 _state[_in_edge] = _state[_in_edge]; … … 1048 1061 } 1049 1062 1050 // Update _thread and _parent node maps1063 // Update _thread and _parent vectors 1051 1064 void updateThreadParent() { 1052 Nodeu;1065 int u; 1053 1066 v_out = _parent[u_out]; 1054 1067 … … 1106 1119 } 1107 1120 1108 // Update _pred _edge and _forward node maps.1121 // Update _pred and _forward vectors 1109 1122 void updatePredEdge() { 1110 Nodeu = u_out, v;1123 int u = u_out, v; 1111 1124 while (u != u_in) { 1112 1125 v = _parent[u]; 1113 _pred _edge[u] = _pred_edge[v];1126 _pred[u] = _pred[v]; 1114 1127 _forward[u] = !_forward[v]; 1115 1128 u = v; 1116 1129 } 1117 _pred _edge[u_in] = _in_edge;1118 _forward[u_in] = (u_in == _ graph.source(_in_edge));1119 } 1120 1121 // Update _depth and _potential node maps1130 _pred[u_in] = _in_edge; 1131 _forward[u_in] = (u_in == _source[_in_edge]); 1132 } 1133 1134 // Update _depth and _potential vectors 1122 1135 void updateDepthPotential() { 1123 1136 _depth[u_in] = _depth[v_in] + 1; 1124 1137 Cost sigma = _forward[u_in] ? 1125 _p otential[v_in]  _potential[u_in]  _cost[_pred_edge[u_in]] :1126 _p otential[v_in]  _potential[u_in] + _cost[_pred_edge[u_in]];1127 _p otential[u_in] += sigma;1128 for( Node u = _thread[u_in]; _parent[u] != INVALID; u = _thread[u]) {1138 _pi[v_in]  _pi[u_in]  _cost[_pred[u_in]] : 1139 _pi[v_in]  _pi[u_in] + _cost[_pred[u_in]]; 1140 _pi[u_in] += sigma; 1141 for(int u = _thread[u_in]; _parent[u] != 1; u = _thread[u]) { 1129 1142 _depth[u] = _depth[_parent[u]] + 1; 1130 1143 if (_depth[u] <= _depth[u_in]) break; 1131 _p otential[u] += sigma;1144 _pi[u] += sigma; 1132 1145 } 1133 1146 } … … 1153 1166 template<class PivotRuleImplementation> 1154 1167 bool start() { 1155 PivotRuleImplementation pivot(*this , _edges);1168 PivotRuleImplementation pivot(*this); 1156 1169 1157 1170 // Execute the network simplex algorithm … … 1159 1172 findJoinNode(); 1160 1173 bool change = findLeavingEdge(); 1161 changeFlow s(change);1174 changeFlow(change); 1162 1175 if (change) { 1163 1176 updateThreadParent(); … … 1168 1181 1169 1182 // Check if the flow amount equals zero on all the artificial edges 1170 for ( InEdgeIt e(_graph, _root); e != INVALID; ++e)1183 for (int e = _edge_num; e != _edge_num + _node_num; ++e) { 1171 1184 if (_flow[e] > 0) return false; 1172 for (OutEdgeIt e(_graph, _root); e != INVALID; ++e) 1173 if (_flow[e] > 0) return false; 1185 } 1174 1186 1175 1187 // Copy flow values to _flow_result 1176 if (_lower) { 1177 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) 1178 (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]]; 1188 if (_orig_lower) { 1189 for (int i = 0; i != _edge_num; ++i) { 1190 Edge e = _edge[i]; 1191 (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i]; 1192 } 1179 1193 } else { 1180 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) 1181 (*_flow_result)[e] = _flow[_edge_ref[e]]; 1194 for (int i = 0; i != _edge_num; ++i) { 1195 (*_flow_result)[_edge[i]] = _flow[i]; 1196 } 1182 1197 } 1183 1198 // Copy potential values to _potential_result 1184 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) 1185 (*_potential_result)[n] = _potential[_node_ref[n]]; 1199 for (int i = 0; i != _node_num; ++i) { 1200 (*_potential_result)[_node[i]] = _pi[i]; 1201 } 1186 1202 1187 1203 return true;
Note: See TracChangeset
for help on using the changeset viewer.