Changeset 2575:e866e288cba6 in lemon0.x for lemon
 Timestamp:
 02/18/08 04:32:06 (12 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@3457
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/network_simplex.h
r2569 r2575 23 23 /// 24 24 /// \file 25 /// \brief The network simplex algorithm for finding a minimum cost flow. 26 25 /// \brief Network simplex algorithm for finding a minimum cost flow. 26 27 #include <vector> 27 28 #include <limits> 29 28 30 #include <lemon/graph_adaptor.h> 29 31 #include <lemon/graph_utils.h> 30 32 #include <lemon/smart_graph.h> 31 32 /// \brief The pivot rule used in the algorithm. 33 //#define FIRST_ELIGIBLE_PIVOT 34 //#define BEST_ELIGIBLE_PIVOT 35 #define EDGE_BLOCK_PIVOT 36 //#define CANDIDATE_LIST_PIVOT 37 //#define SORTED_LIST_PIVOT 38 39 //#define _DEBUG_ITER_ 40 41 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 48 49 #ifdef EDGE_BLOCK_PIVOT 50 #include <lemon/math.h> 51 #define MIN_BLOCK_SIZE 10 52 #endif 53 54 #ifdef CANDIDATE_LIST_PIVOT 55 #include <vector> 56 #define LIST_LENGTH_DIV 50 57 #define MINOR_LIMIT_DIV 200 58 #endif 59 60 #ifdef SORTED_LIST_PIVOT 61 #include <vector> 62 #include <algorithm> 63 #define LIST_LENGTH_DIV 100 64 #define LOWER_DIV 4 65 #endif 33 #include <lemon/math.h> 66 34 67 35 namespace lemon { … … 76 44 /// finding a minimum cost flow. 77 45 /// 78 /// \ param Graph The directed graph type the algorithm runs on.79 /// \ param LowerMap The type of the lower bound map.80 /// \ param CapacityMap The type of the capacity (upper bound) map.81 /// \ param CostMap The type of the cost (length) map.82 /// \ param SupplyMap The type of the supply map.46 /// \tparam Graph The directed graph type the algorithm runs on. 47 /// \tparam LowerMap The type of the lower bound map. 48 /// \tparam CapacityMap The type of the capacity (upper bound) map. 49 /// \tparam CostMap The type of the cost (length) map. 50 /// \tparam SupplyMap The type of the supply map. 83 51 /// 84 52 /// \warning 85 ///  Edge capacities and costs should be nonnegative integers. 86 /// However \c CostMap::Value should be signed type. 87 ///  Supply values should be signed integers. 88 ///  \c LowerMap::Value must be convertible to 89 /// \c CapacityMap::Value and \c CapacityMap::Value must be 90 /// convertible to \c SupplyMap::Value. 53 ///  Edge capacities and costs should be \e nonnegative \e integers. 54 ///  Supply values should be \e signed \e integers. 55 ///  \c LowerMap::Value must be convertible to \c CapacityMap::Value. 56 ///  \c CapacityMap::Value and \c SupplyMap::Value must be 57 /// convertible to each other. 58 ///  All value types must be convertible to \c CostMap::Value, which 59 /// must be signed type. 60 /// 61 /// \note \ref NetworkSimplex provides six different pivot rule 62 /// implementations that significantly affect the efficiency of the 63 /// algorithm. 64 /// By default a combined pivot rule is used, which is the fastest 65 /// implementation according to our benchmark tests. 66 /// Another pivot rule can be selected using \ref run() function 67 /// with the proper parameter. 91 68 /// 92 69 /// \author Peter Kovacs … … 94 71 template < typename Graph, 95 72 typename LowerMap = typename Graph::template EdgeMap<int>, 96 typename CapacityMap = LowerMap,73 typename CapacityMap = typename Graph::template EdgeMap<int>, 97 74 typename CostMap = typename Graph::template EdgeMap<int>, 98 typename SupplyMap = typename Graph::template NodeMap 99 <typename CapacityMap::Value> > 75 typename SupplyMap = typename Graph::template NodeMap<int> > 100 76 class NetworkSimplex 101 77 { 102 typedef typename LowerMap::Value Lower;103 78 typedef typename CapacityMap::Value Capacity; 104 79 typedef typename CostMap::Value Cost; … … 108 83 GRAPH_TYPEDEFS(typename SGraph); 109 84 110 typedef typename SGraph::template EdgeMap<Lower> SLowerMap;111 85 typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap; 112 86 typedef typename SGraph::template EdgeMap<Cost> SCostMap; … … 130 104 typedef typename Graph::template NodeMap<Cost> PotentialMap; 131 105 132 protected: 133 106 public: 107 108 /// Enum type to select the pivot rule used by \ref run(). 109 enum PivotRuleEnum { 110 FIRST_ELIGIBLE_PIVOT, 111 BEST_ELIGIBLE_PIVOT, 112 BLOCK_SEARCH_PIVOT, 113 LIMITED_SEARCH_PIVOT, 114 CANDIDATE_LIST_PIVOT, 115 COMBINED_PIVOT 116 }; 117 118 private: 119 120 /// \brief Map adaptor class for handling reduced edge costs. 121 /// 134 122 /// Map adaptor class for handling reduced edge costs. 135 123 class ReducedCostMap : public MapBase<Edge, Cost> … … 137 125 private: 138 126 139 const SGraph & gr;140 const SCostMap & cost_map;141 const SPotentialMap & pot_map;127 const SGraph &_gr; 128 const SCostMap &_cost_map; 129 const SPotentialMap &_pot_map; 142 130 143 131 public: 144 132 145 ReducedCostMap( const SGraph &_gr, 146 const SCostMap &_cm, 147 const SPotentialMap &_pm ) : 148 gr(_gr), cost_map(_cm), pot_map(_pm) {} 149 133 ///\e 134 ReducedCostMap( const SGraph &gr, 135 const SCostMap &cost_map, 136 const SPotentialMap &pot_map ) : 137 _gr(gr), _cost_map(cost_map), _pot_map(pm) {} 138 139 ///\e 150 140 Cost operator[](const Edge &e) const { 151 return cost_map[e]  pot_map[gr.source(e)] + pot_map[gr.target(e)]; 141 return _cost_map[e] + _pot_map[_gr.source(e)] 142  _pot_map[_gr.target(e)]; 152 143 } 153 144 154 145 }; //class ReducedCostMap 155 146 156 protected: 157 158 /// The directed graph the algorithm runs on. 159 SGraph graph; 160 /// The original graph. 161 const Graph &graph_ref; 162 /// The original lower bound map. 163 const LowerMap *lower; 164 /// The capacity map. 165 SCapacityMap capacity; 166 /// The cost map. 167 SCostMap cost; 168 /// The supply map. 169 SSupplyMap supply; 170 /// The reduced cost map. 171 ReducedCostMap red_cost; 172 bool valid_supply; 173 174 /// The edge map of the current flow. 175 SCapacityMap flow; 176 /// The potential node map. 177 SPotentialMap potential; 178 FlowMap flow_result; 179 PotentialMap potential_result; 180 181 /// Node reference for the original graph. 182 NodeRefMap node_ref; 183 /// Edge reference for the original graph. 184 EdgeRefMap edge_ref; 185 186 /// The \c depth node map of the spanning tree structure. 187 IntNodeMap depth; 188 /// The \c parent node map of the spanning tree structure. 189 NodeNodeMap parent; 190 /// The \c pred_edge node map of the spanning tree structure. 191 EdgeNodeMap pred_edge; 192 /// The \c thread node map of the spanning tree structure. 193 NodeNodeMap thread; 194 /// The \c forward node map of the spanning tree structure. 195 BoolNodeMap forward; 196 /// The state edge map. 197 IntEdgeMap state; 198 /// The root node of the starting spanning tree. 199 Node root; 147 private: 148 149 /// \brief Implementation of the "First Eligible" pivot rule for the 150 /// \ref NetworkSimplex "network simplex" algorithm. 151 /// 152 /// This class implements the "First Eligible" pivot rule 153 /// for the \ref NetworkSimplex "network simplex" algorithm. 154 class FirstEligiblePivotRule 155 { 156 private: 157 158 NetworkSimplex &_ns; 159 EdgeIt _next_edge; 160 161 public: 162 163 /// Constructor. 164 FirstEligiblePivotRule(NetworkSimplex &ns) : 165 _ns(ns), _next_edge(ns._graph) {} 166 167 /// Finds the next entering edge. 168 bool findEnteringEdge() { 169 for (EdgeIt e = _next_edge; e != INVALID; ++e) { 170 if (_ns._state[e] * _ns._red_cost[e] < 0) { 171 _ns._in_edge = e; 172 _next_edge = ++e; 173 return true; 174 } 175 } 176 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { 177 if (_ns._state[e] * _ns._red_cost[e] < 0) { 178 _ns._in_edge = e; 179 _next_edge = ++e; 180 return true; 181 } 182 } 183 return false; 184 } 185 }; //class FirstEligiblePivotRule 186 187 /// \brief Implementation of the "Best Eligible" pivot rule for the 188 /// \ref NetworkSimplex "network simplex" algorithm. 189 /// 190 /// This class implements the "Best Eligible" pivot rule 191 /// for the \ref NetworkSimplex "network simplex" algorithm. 192 class BestEligiblePivotRule 193 { 194 private: 195 196 NetworkSimplex &_ns; 197 198 public: 199 200 /// Constructor. 201 BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {} 202 203 /// Finds the next entering edge. 204 bool findEnteringEdge() { 205 Cost min = 0; 206 for (EdgeIt e(_ns._graph); e != INVALID; ++e) { 207 if (_ns._state[e] * _ns._red_cost[e] < min) { 208 min = _ns._state[e] * _ns._red_cost[e]; 209 _ns._in_edge = e; 210 } 211 } 212 return min < 0; 213 } 214 }; //class BestEligiblePivotRule 215 216 /// \brief Implementation of the "Block Search" pivot rule for the 217 /// \ref NetworkSimplex "network simplex" algorithm. 218 /// 219 /// This class implements the "Block Search" pivot rule 220 /// for the \ref NetworkSimplex "network simplex" algorithm. 221 class BlockSearchPivotRule 222 { 223 private: 224 225 NetworkSimplex &_ns; 226 EdgeIt _next_edge, _min_edge; 227 int _block_size; 228 229 static const int MIN_BLOCK_SIZE = 10; 230 231 public: 232 233 /// Constructor. 234 BlockSearchPivotRule(NetworkSimplex &ns) : 235 _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) 236 { 237 _block_size = 2 * int(sqrt(countEdges(_ns._graph))); 238 if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE; 239 } 240 241 /// Finds the next entering edge. 242 bool findEnteringEdge() { 243 Cost curr, min = 0; 244 int cnt = 0; 245 for (EdgeIt e = _next_edge; e != INVALID; ++e) { 246 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { 247 min = curr; 248 _min_edge = e; 249 } 250 if (++cnt == _block_size) { 251 if (min < 0) break; 252 cnt = 0; 253 } 254 } 255 if (min == 0) { 256 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { 257 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { 258 min = curr; 259 _min_edge = e; 260 } 261 if (++cnt == _block_size) { 262 if (min < 0) break; 263 cnt = 0; 264 } 265 } 266 } 267 _ns._in_edge = _min_edge; 268 _next_edge = ++_min_edge; 269 return min < 0; 270 } 271 }; //class BlockSearchPivotRule 272 273 /// \brief Implementation of the "Limited Search" pivot rule for the 274 /// \ref NetworkSimplex "network simplex" algorithm. 275 /// 276 /// This class implements the "Limited Search" pivot rule 277 /// for the \ref NetworkSimplex "network simplex" algorithm. 278 class LimitedSearchPivotRule 279 { 280 private: 281 282 NetworkSimplex &_ns; 283 EdgeIt _next_edge, _min_edge; 284 int _sample_size; 285 286 static const int MIN_SAMPLE_SIZE = 10; 287 static const double SAMPLE_SIZE_FACTOR = 0.0015; 288 289 public: 290 291 /// Constructor. 292 LimitedSearchPivotRule(NetworkSimplex &ns) : 293 _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) 294 { 295 _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph)); 296 if (_sample_size < MIN_SAMPLE_SIZE) 297 _sample_size = MIN_SAMPLE_SIZE; 298 } 299 300 /// Finds the next entering edge. 301 bool findEnteringEdge() { 302 Cost curr, min = 0; 303 int cnt = 0; 304 for (EdgeIt e = _next_edge; e != INVALID; ++e) { 305 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { 306 min = curr; 307 _min_edge = e; 308 } 309 if (curr < 0 && ++cnt == _sample_size) break; 310 } 311 if (min == 0) { 312 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { 313 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { 314 min = curr; 315 _min_edge = e; 316 } 317 if (curr < 0 && ++cnt == _sample_size) break; 318 } 319 } 320 _ns._in_edge = _min_edge; 321 _next_edge = ++_min_edge; 322 return min < 0; 323 } 324 }; //class LimitedSearchPivotRule 325 326 /// \brief Implementation of the "Candidate List" pivot rule for the 327 /// \ref NetworkSimplex "network simplex" algorithm. 328 /// 329 /// This class implements the "Candidate List" pivot rule 330 /// for the \ref NetworkSimplex "network simplex" algorithm. 331 class CandidateListPivotRule 332 { 333 private: 334 335 NetworkSimplex &_ns; 336 337 // The list of candidate edges. 338 std::vector<Edge> _candidates; 339 // The maximum length of the edge list. 340 int _list_length; 341 // The maximum number of minor iterations between two major 342 // itarations. 343 int _minor_limit; 344 345 int _minor_count; 346 EdgeIt _next_edge; 347 348 static const double LIST_LENGTH_FACTOR = 0.002; 349 static const double MINOR_LIMIT_FACTOR = 0.1; 350 static const int MIN_LIST_LENGTH = 10; 351 static const int MIN_MINOR_LIMIT = 2; 352 353 public: 354 355 /// Constructor. 356 CandidateListPivotRule(NetworkSimplex &ns) : 357 _ns(ns), _next_edge(ns._graph) 358 { 359 int edge_num = countEdges(_ns._graph); 360 _minor_count = 0; 361 _list_length = int(edge_num * LIST_LENGTH_FACTOR); 362 if (_list_length < MIN_LIST_LENGTH) 363 _list_length = MIN_LIST_LENGTH; 364 _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR); 365 if (_minor_limit < MIN_MINOR_LIMIT) 366 _minor_limit = MIN_MINOR_LIMIT; 367 } 368 369 /// Finds the next entering edge. 370 bool findEnteringEdge() { 371 Cost min, curr; 372 if (_minor_count < _minor_limit && _candidates.size() > 0) { 373 // Minor iteration 374 ++_minor_count; 375 Edge e; 376 min = 0; 377 for (int i = 0; i < int(_candidates.size()); ++i) { 378 e = _candidates[i]; 379 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { 380 min = curr; 381 _ns._in_edge = e; 382 } 383 } 384 if (min < 0) return true; 385 } 386 387 // Major iteration 388 _candidates.clear(); 389 EdgeIt e = _next_edge; 390 min = 0; 391 for ( ; e != INVALID; ++e) { 392 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { 393 _candidates.push_back(e); 394 if (curr < min) { 395 min = curr; 396 _ns._in_edge = e; 397 } 398 if (int(_candidates.size()) == _list_length) break; 399 } 400 } 401 if (int(_candidates.size()) < _list_length) { 402 for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) { 403 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { 404 _candidates.push_back(e); 405 if (curr < min) { 406 min = curr; 407 _ns._in_edge = e; 408 } 409 if (int(_candidates.size()) == _list_length) break; 410 } 411 } 412 } 413 if (_candidates.size() == 0) return false; 414 _minor_count = 1; 415 _next_edge = ++e; 416 return true; 417 } 418 }; //class CandidateListPivotRule 419 420 private: 421 422 // State constant for edges at their lower bounds. 423 static const int STATE_LOWER = 1; 424 // State constant for edges in the spanning tree. 425 static const int STATE_TREE = 0; 426 // State constant for edges at their upper bounds. 427 static const int STATE_UPPER = 1; 428 429 // Constant for the combined pivot rule. 430 static const int COMBINED_PIVOT_MAX_DEG = 5; 431 432 private: 433 434 // The directed graph the algorithm runs on 435 SGraph _graph; 436 // The original graph 437 const Graph &_graph_ref; 438 // The original lower bound map 439 const LowerMap *_lower; 440 // The capacity map 441 SCapacityMap _capacity; 442 // The cost map 443 SCostMap _cost; 444 // The supply map 445 SSupplyMap _supply; 446 bool _valid_supply; 447 448 // Edge map of the current flow 449 SCapacityMap _flow; 450 // Node map of the current potentials 451 SPotentialMap _potential; 452 453 // The depth node map of the spanning tree structure 454 IntNodeMap _depth; 455 // The parent node map of the spanning tree structure 456 NodeNodeMap _parent; 457 // The pred_edge node map of the spanning tree structure 458 EdgeNodeMap _pred_edge; 459 // The thread node map of the spanning tree structure 460 NodeNodeMap _thread; 461 // The forward node map of the spanning tree structure 462 BoolNodeMap _forward; 463 // The state edge map 464 IntEdgeMap _state; 465 // The root node of the starting spanning tree 466 Node _root; 467 468 // The reduced cost map 469 ReducedCostMap _red_cost; 470 471 // Members for handling the original graph 472 FlowMap _flow_result; 473 PotentialMap _potential_result; 474 NodeRefMap _node_ref; 475 EdgeRefMap _edge_ref; 200 476 201 477 // The entering edge of the current pivot iteration. 202 Edge in_edge; 478 Edge _in_edge; 479 203 480 // Temporary nodes used in the current pivot iteration. 204 481 Node join, u_in, v_in, u_out, v_out; … … 209 486 Capacity delta; 210 487 211 #ifdef EDGE_BLOCK_PIVOT212 /// The size of blocks for the "Edge Block" pivot rule.213 int block_size;214 #endif215 #ifdef CANDIDATE_LIST_PIVOT216 /// \brief The list of candidate edges for the "Candidate List"217 /// pivot rule.218 std::vector<Edge> candidates;219 /// \brief The maximum length of the edge list for the220 /// "Candidate List" pivot rule.221 int list_length;222 /// \brief The maximum number of minor iterations between two major223 /// itarations.224 int minor_limit;225 /// \brief The number of minor iterations.226 int minor_count;227 #endif228 #ifdef SORTED_LIST_PIVOT229 /// \brief The list of candidate edges for the "Sorted List"230 /// pivot rule.231 std::vector<Edge> candidates;232 /// \brief The maximum length of the edge list for the233 /// "Sorted List" pivot rule.234 int list_length;235 int list_index;236 #endif237 238 488 public : 239 489 … … 242 492 /// General constructor of the class (with lower bounds). 243 493 /// 244 /// \param _graph The directed graph the algorithm runs on. 245 /// \param _lower The lower bounds of the edges. 246 /// \param _capacity The capacities (upper bounds) of the edges. 247 /// \param _cost The cost (length) values of the edges. 248 /// \param _supply The supply values of the nodes (signed). 249 NetworkSimplex( const Graph &_graph, 250 const LowerMap &_lower, 251 const CapacityMap &_capacity, 252 const CostMap &_cost, 253 const SupplyMap &_supply ) : 254 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), 255 supply(graph), flow(graph), flow_result(_graph), potential(graph), 256 potential_result(_graph), depth(graph), parent(graph), 257 pred_edge(graph), thread(graph), forward(graph), state(graph), 258 node_ref(graph_ref), edge_ref(graph_ref), 259 red_cost(graph, cost, potential) 494 /// \param graph The directed graph the algorithm runs on. 495 /// \param lower The lower bounds of the edges. 496 /// \param capacity The capacities (upper bounds) of the edges. 497 /// \param cost The cost (length) values of the edges. 498 /// \param supply The supply values of the nodes (signed). 499 NetworkSimplex( const Graph &graph, 500 const LowerMap &lower, 501 const CapacityMap &capacity, 502 const CostMap &cost, 503 const SupplyMap &supply ) : 504 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), 505 _cost(_graph), _supply(_graph), _flow(_graph), 506 _potential(_graph), _depth(_graph), _parent(_graph), 507 _pred_edge(_graph), _thread(_graph), _forward(_graph), 508 _state(_graph), _red_cost(_graph, _cost, _potential), 509 _flow_result(graph), _potential_result(graph), 510 _node_ref(graph), _edge_ref(graph) 260 511 { 261 512 // Checking the sum of supply values 262 513 Supply sum = 0; 263 for (typename Graph::NodeIt n( graph_ref); n != INVALID; ++n)264 sum += _supply[n];265 if (!( valid_supply = sum == 0)) return;266 267 // Copying graph_ref tograph268 graph.reserveNode(countNodes(graph_ref) + 1);269 graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));270 copyGraph( graph,graph_ref)271 .edgeMap( cost, _cost)272 .nodeRef( node_ref)273 .edgeRef( edge_ref)514 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) 515 sum += supply[n]; 516 if (!(_valid_supply = sum == 0)) return; 517 518 // Copying _graph_ref to _graph 519 _graph.reserveNode(countNodes(_graph_ref) + 1); 520 _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); 521 copyGraph(_graph, _graph_ref) 522 .edgeMap(_cost, cost) 523 .nodeRef(_node_ref) 524 .edgeRef(_edge_ref) 274 525 .run(); 275 526 276 527 // Removing nonzero lower bounds 277 for (typename Graph::EdgeIt e( graph_ref); e != INVALID; ++e) {278 capacity[edge_ref[e]] = _capacity[e]  _lower[e];279 } 280 for (typename Graph::NodeIt n( graph_ref); n != INVALID; ++n) {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;528 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { 529 _capacity[_edge_ref[e]] = capacity[e]  lower[e]; 530 } 531 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { 532 Supply s = supply[n]; 533 for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) 534 s += lower[e]; 535 for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) 536 s = lower[e]; 537 _supply[_node_ref[n]] = s; 287 538 } 288 539 } … … 292 543 /// General constructor of the class (without lower bounds). 293 544 /// 294 /// \param _graph The directed graph the algorithm runs on. 295 /// \param _capacity The capacities (upper bounds) of the edges. 296 /// \param _cost The cost (length) values of the edges. 297 /// \param _supply The supply values of the nodes (signed). 298 NetworkSimplex( const Graph &_graph, 299 const CapacityMap &_capacity, 300 const CostMap &_cost, 301 const SupplyMap &_supply ) : 302 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), 303 supply(graph), flow(graph), flow_result(_graph), potential(graph), 304 potential_result(_graph), depth(graph), parent(graph), 305 pred_edge(graph), thread(graph), forward(graph), state(graph), 306 node_ref(graph_ref), edge_ref(graph_ref), 307 red_cost(graph, cost, potential) 545 /// \param graph The directed graph the algorithm runs on. 546 /// \param capacity The capacities (upper bounds) of the edges. 547 /// \param cost The cost (length) values of the edges. 548 /// \param supply The supply values of the nodes (signed). 549 NetworkSimplex( const Graph &graph, 550 const CapacityMap &capacity, 551 const CostMap &cost, 552 const SupplyMap &supply ) : 553 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), 554 _cost(_graph), _supply(_graph), _flow(_graph), 555 _potential(_graph), _depth(_graph), _parent(_graph), 556 _pred_edge(_graph), _thread(_graph), _forward(_graph), 557 _state(_graph), _red_cost(_graph, _cost, _potential), 558 _flow_result(graph), _potential_result(graph), 559 _node_ref(graph), _edge_ref(graph) 308 560 { 309 561 // Checking the sum of supply values 310 562 Supply sum = 0; 311 for (typename Graph::NodeIt n( graph_ref); n != INVALID; ++n)312 sum += _supply[n];313 if (!( valid_supply = sum == 0)) return;314 315 // Copying graph_ref to graph316 copyGraph( graph,graph_ref)317 .edgeMap( capacity, _capacity)318 .edgeMap( cost, _cost)319 .nodeMap( supply, _supply)320 .nodeRef( node_ref)321 .edgeRef( edge_ref)563 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) 564 sum += supply[n]; 565 if (!(_valid_supply = sum == 0)) return; 566 567 // Copying _graph_ref to graph 568 copyGraph(_graph, _graph_ref) 569 .edgeMap(_capacity, capacity) 570 .edgeMap(_cost, cost) 571 .nodeMap(_supply, supply) 572 .nodeRef(_node_ref) 573 .edgeRef(_edge_ref) 322 574 .run(); 323 575 } … … 327 579 /// Simple constructor of the class (with lower bounds). 328 580 /// 329 /// \param _graph The directed graph the algorithm runs on. 330 /// \param _lower The lower bounds of the edges. 331 /// \param _capacity The capacities (upper bounds) of the edges. 332 /// \param _cost The cost (length) values of the edges. 333 /// \param _s The source node. 334 /// \param _t The target node. 335 /// \param _flow_value The required amount of flow from node \c _s 336 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). 337 NetworkSimplex( const Graph &_graph, 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 ) : 344 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), 345 supply(graph), flow(graph), flow_result(_graph), potential(graph), 346 potential_result(_graph), depth(graph), parent(graph), 347 pred_edge(graph), thread(graph), forward(graph), state(graph), 348 node_ref(graph_ref), edge_ref(graph_ref), 349 red_cost(graph, cost, potential) 581 /// \param graph The directed graph the algorithm runs on. 582 /// \param lower The lower bounds of the edges. 583 /// \param capacity The capacities (upper bounds) of the edges. 584 /// \param cost The cost (length) values of the edges. 585 /// \param s The source node. 586 /// \param t The target node. 587 /// \param flow_value The required amount of flow from node \c s 588 /// to node \c t (i.e. the supply of \c s and the demand of \c t). 589 NetworkSimplex( const Graph &graph, 590 const LowerMap &lower, 591 const CapacityMap &capacity, 592 const CostMap &cost, 593 typename Graph::Node s, 594 typename Graph::Node t, 595 typename SupplyMap::Value flow_value ) : 596 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), 597 _cost(_graph), _supply(_graph), _flow(_graph), 598 _potential(_graph), _depth(_graph), _parent(_graph), 599 _pred_edge(_graph), _thread(_graph), _forward(_graph), 600 _state(_graph), _red_cost(_graph, _cost, _potential), 601 _flow_result(graph), _potential_result(graph), 602 _node_ref(graph), _edge_ref(graph) 350 603 { 351 // Copying graph_ref to graph352 copyGraph( graph,graph_ref)353 .edgeMap( cost, _cost)354 .nodeRef( node_ref)355 .edgeRef( edge_ref)604 // Copying _graph_ref to graph 605 copyGraph(_graph, _graph_ref) 606 .edgeMap(_cost, cost) 607 .nodeRef(_node_ref) 608 .edgeRef(_edge_ref) 356 609 .run(); 357 610 358 611 // Removing nonzero lower bounds 359 for (typename Graph::EdgeIt e( graph_ref); e != INVALID; ++e) {360 capacity[edge_ref[e]] = _capacity[e]  _lower[e];361 } 362 for (typename Graph::NodeIt n( graph_ref); n != INVALID; ++n) {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;371 } 372 valid_supply = true;612 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { 613 _capacity[_edge_ref[e]] = capacity[e]  lower[e]; 614 } 615 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { 616 Supply sum = 0; 617 if (n == s) sum = flow_value; 618 if (n == t) sum = flow_value; 619 for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) 620 sum += lower[e]; 621 for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) 622 sum = lower[e]; 623 _supply[_node_ref[n]] = sum; 624 } 625 _valid_supply = true; 373 626 } 374 627 … … 377 630 /// Simple constructor of the class (without lower bounds). 378 631 /// 379 /// \param _graph The directed graph the algorithm runs on. 380 /// \param _capacity The capacities (upper bounds) of the edges. 381 /// \param _cost The cost (length) values of the edges. 382 /// \param _s The source node. 383 /// \param _t The target node. 384 /// \param _flow_value The required amount of flow from node \c _s 385 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). 386 NetworkSimplex( const Graph &_graph, 387 const CapacityMap &_capacity, 388 const CostMap &_cost, 389 typename Graph::Node _s, 390 typename Graph::Node _t, 391 typename SupplyMap::Value _flow_value ) : 392 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), 393 supply(graph, 0), flow(graph), flow_result(_graph), potential(graph), 394 potential_result(_graph), depth(graph), parent(graph), 395 pred_edge(graph), thread(graph), forward(graph), state(graph), 396 node_ref(graph_ref), edge_ref(graph_ref), 397 red_cost(graph, cost, potential) 632 /// \param graph The directed graph the algorithm runs on. 633 /// \param capacity The capacities (upper bounds) of the edges. 634 /// \param cost The cost (length) values of the edges. 635 /// \param s The source node. 636 /// \param t The target node. 637 /// \param flow_value The required amount of flow from node \c s 638 /// to node \c t (i.e. the supply of \c s and the demand of \c t). 639 NetworkSimplex( const Graph &graph, 640 const CapacityMap &capacity, 641 const CostMap &cost, 642 typename Graph::Node s, 643 typename Graph::Node t, 644 typename SupplyMap::Value flow_value ) : 645 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), 646 _cost(_graph), _supply(_graph, 0), _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), 650 _flow_result(graph), _potential_result(graph), 651 _node_ref(graph), _edge_ref(graph) 398 652 { 399 // Copying graph_ref to graph400 copyGraph( graph,graph_ref)401 .edgeMap( capacity, _capacity)402 .edgeMap( cost, _cost)403 .nodeRef( node_ref)404 .edgeRef( edge_ref)653 // Copying _graph_ref to graph 654 copyGraph(_graph, _graph_ref) 655 .edgeMap(_capacity, capacity) 656 .edgeMap(_cost, cost) 657 .nodeRef(_node_ref) 658 .edgeRef(_edge_ref) 405 659 .run(); 406 supply[node_ref[_s]] = _flow_value;407 supply[node_ref[_t]] = _flow_value;408 valid_supply = true;660 _supply[_node_ref[s]] = flow_value; 661 _supply[_node_ref[t]] = flow_value; 662 _valid_supply = true; 409 663 } 410 664 … … 413 667 /// Runs the algorithm. 414 668 /// 669 /// \param pivot_rule The pivot rule that is used during the 670 /// algorithm. 671 /// 672 /// The available pivot rules: 673 /// 674 ///  FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in 675 /// a wraparound fashion in every iteration 676 /// (\ref FirstEligiblePivotRule). 677 /// 678 ///  BEST_ELIGIBLE_PIVOT The best eligible edge is selected in 679 /// every iteration (\ref BestEligiblePivotRule). 680 /// 681 ///  BLOCK_SEARCH_PIVOT A specified number of edges are examined in 682 /// every iteration in a wraparound fashion and the best eligible 683 /// edge is selected from this block (\ref BlockSearchPivotRule). 684 /// 685 ///  LIMITED_SEARCH_PIVOT A specified number of eligible edges are 686 /// examined in every iteration in a wraparound fashion and the best 687 /// one is selected from them (\ref LimitedSearchPivotRule). 688 /// 689 ///  CANDIDATE_LIST_PIVOT In major iterations a candidate list is 690 /// built from eligible edges and it is used for edge selection in 691 /// the following minor iterations (\ref CandidateListPivotRule). 692 /// 693 ///  COMBINED_PIVOT This is a combined version of the two fastest 694 /// pivot rules. 695 /// For rather sparse graphs \ref LimitedSearchPivotRule 696 /// "Limited Search" implementation is used, otherwise 697 /// \ref BlockSearchPivotRule "Block Search" pivot rule is used. 698 /// According to our benchmark tests this combined method is the 699 /// most efficient. 700 /// 415 701 /// \return \c true if a feasible flow can be found. 416 bool run() { 417 return init() && start(); 418 } 419 420 /// \brief Returns a const reference to the flow map. 421 /// 422 /// Returns a const reference to the flow map. 702 bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) { 703 return init() && start(pivot_rule); 704 } 705 706 /// \brief Returns a const reference to the edge map storing the 707 /// found flow. 708 /// 709 /// Returns a const reference to the edge map storing the found flow. 423 710 /// 424 711 /// \pre \ref run() must be called before using this function. 425 712 const FlowMap& flowMap() const { 426 return flow_result;427 } 428 429 /// \brief Returns a const reference to the potential map (the dual430 /// solution).431 /// 432 /// Returns a const reference to the potential map (the dual433 /// solution).713 return _flow_result; 714 } 715 716 /// \brief Returns a const reference to the node map storing the 717 /// found potentials (the dual solution). 718 /// 719 /// Returns a const reference to the node map storing the found 720 /// potentials (the dual solution). 434 721 /// 435 722 /// \pre \ref run() must be called before using this function. 436 723 const PotentialMap& potentialMap() const { 437 return potential_result;724 return _potential_result; 438 725 } 439 726 … … 446 733 Cost totalCost() const { 447 734 Cost c = 0; 448 for (typename Graph::EdgeIt e( graph_ref); e != INVALID; ++e)449 c += flow_result[e] * cost[edge_ref[e]];735 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) 736 c += _flow_result[e] * _cost[_edge_ref[e]]; 450 737 return c; 451 738 } 452 739 453 pr otected:740 private: 454 741 455 742 /// \brief Extends the underlaying graph and initializes all the 456 743 /// node and edge maps. 457 744 bool init() { 458 if (! valid_supply) return false;745 if (!_valid_supply) return false; 459 746 460 747 // Initializing state and flow maps 461 for (EdgeIt e( graph); e != INVALID; ++e) {462 flow[e] = 0;463 state[e] =LOWER;748 for (EdgeIt e(_graph); e != INVALID; ++e) { 749 _flow[e] = 0; 750 _state[e] = STATE_LOWER; 464 751 } 465 752 466 753 // Adding an artificial root node to the graph 467 root =graph.addNode();468 parent[root] = INVALID;469 pred_edge[root] = INVALID;470 depth[root] = 0;471 supply[root] = 0;472 potential[root] = 0;754 _root = _graph.addNode(); 755 _parent[_root] = INVALID; 756 _pred_edge[_root] = INVALID; 757 _depth[_root] = 0; 758 _supply[_root] = 0; 759 _potential[_root] = 0; 473 760 474 761 // Adding artificial edges to the graph and initializing the node 475 762 // maps of the spanning tree data structure 476 Supply sum = 0; 477 Node last = root; 763 Node last = _root; 478 764 Edge e; 479 765 Cost max_cost = std::numeric_limits<Cost>::max() / 4; 480 for (NodeIt u( graph); u != INVALID; ++u) {481 if (u == root) continue;482 thread[last] = u;766 for (NodeIt u(_graph); u != INVALID; ++u) { 767 if (u == _root) continue; 768 _thread[last] = u; 483 769 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; 770 _parent[u] = _root; 771 _depth[u] = 1; 772 if (_supply[u] >= 0) { 773 e = _graph.addEdge(u, _root); 774 _flow[e] = _supply[u]; 775 _forward[u] = true; 776 _potential[u] = max_cost; 492 777 } 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; 502 } 503 thread[last] = root; 504 505 #ifdef EDGE_BLOCK_PIVOT 506 // Initializing block_size for the edge block pivot rule 507 int edge_num = countEdges(graph); 508 block_size = 2 * int(sqrt(countEdges(graph))); 509 if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE; 510 #endif 511 #ifdef CANDIDATE_LIST_PIVOT 512 int edge_num = countEdges(graph); 513 minor_count = 0; 514 list_length = edge_num / LIST_LENGTH_DIV; 515 minor_limit = edge_num / MINOR_LIMIT_DIV; 516 #endif 517 #ifdef SORTED_LIST_PIVOT 518 int edge_num = countEdges(graph); 519 list_index = 0; 520 list_length = edge_num / LIST_LENGTH_DIV; 521 #endif 522 523 return sum == 0; 524 } 525 526 #ifdef FIRST_ELIGIBLE_PIVOT 527 /// \brief Finds entering edge according to the "First Eligible" 528 /// pivot rule. 529 bool findEnteringEdge(EdgeIt &next_edge) { 530 for (EdgeIt e = next_edge; e != INVALID; ++e) { 531 if (state[e] * red_cost[e] < 0) { 532 in_edge = e; 533 next_edge = ++e; 534 return true; 535 } 536 } 537 for (EdgeIt e(graph); e != next_edge; ++e) { 538 if (state[e] * red_cost[e] < 0) { 539 in_edge = e; 540 next_edge = ++e; 541 return true; 542 } 543 } 544 return false; 545 } 546 #endif 547 548 #ifdef BEST_ELIGIBLE_PIVOT 549 /// \brief Finds entering edge according to the "Best Eligible" 550 /// pivot rule. 551 bool findEnteringEdge() { 552 Cost min = 0; 553 for (EdgeIt e(graph); e != INVALID; ++e) { 554 if (state[e] * red_cost[e] < min) { 555 min = state[e] * red_cost[e]; 556 in_edge = e; 557 } 558 } 559 return min < 0; 560 } 561 #endif 562 563 #ifdef EDGE_BLOCK_PIVOT 564 /// \brief Finds entering edge according to the "Edge Block" 565 /// pivot rule. 566 bool findEnteringEdge(EdgeIt &next_edge) { 567 // Performing edge block selection 568 Cost curr, min = 0; 569 EdgeIt min_edge(graph); 570 int cnt = 0; 571 for (EdgeIt e = next_edge; e != INVALID; ++e) { 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 } 580 } 581 if (!(min < 0)) { 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 } 592 } 593 in_edge = min_edge; 594 if ((next_edge = ++min_edge) == INVALID) 595 next_edge = EdgeIt(graph); 596 return min < 0; 597 } 598 #endif 599 600 #ifdef CANDIDATE_LIST_PIVOT 601 /// \brief Finds entering edge according to the "Candidate List" 602 /// pivot rule. 603 bool findEnteringEdge() { 604 typedef typename std::vector<Edge>::iterator ListIt; 605 606 if (minor_count >= minor_limit  candidates.size() == 0) { 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; 616 } 617 618 // Minor iteration 619 ++minor_count; 620 Cost min = 0; 621 Edge e; 622 for (int i = 0; i < candidates.size(); ++i) { 623 e = candidates[i]; 624 if (state[e] * red_cost[e] < min) { 625 min = state[e] * red_cost[e]; 626 in_edge = e; 627 } 628 } 778 e = _graph.addEdge(_root, u); 779 _flow[e] = _supply[u]; 780 _forward[u] = false; 781 _potential[u] = max_cost; 782 } 783 _cost[e] = max_cost; 784 _capacity[e] = std::numeric_limits<Capacity>::max(); 785 _state[e] = STATE_TREE; 786 _pred_edge[u] = e; 787 } 788 _thread[last] = _root; 789 629 790 return true; 630 791 } 631 #endif 632 633 #ifdef SORTED_LIST_PIVOT 634 /// \brief Functor class to compare edges during sort of the 635 /// candidate list. 636 class SortFunc 637 { 638 private: 639 const IntEdgeMap &st; 640 const ReducedCostMap &rc; 641 public: 642 SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) : 643 st(_st), rc(_rc) {} 644 bool operator()(const Edge &e1, const Edge &e2) { 645 return st[e1] * rc[e1] < st[e2] * rc[e2]; 646 } 647 }; 648 649 /// \brief Finds entering edge according to the "Sorted List" 650 /// pivot rule. 651 bool findEnteringEdge() { 652 static SortFunc sort_func(state, red_cost); 653 654 // Minor iteration 655 while (list_index < candidates.size()) { 656 in_edge = candidates[list_index++]; 657 if (state[in_edge] * red_cost[in_edge] < 0) return true; 658 } 659 660 // Major iteration 661 candidates.clear(); 662 Cost curr, min = 0; 663 for (EdgeIt e(graph); e != INVALID; ++e) { 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 } 669 } 670 if (candidates.size() == 0) return false; 671 sort(candidates.begin(), candidates.end(), sort_func); 672 in_edge = candidates[0]; 673 list_index = 1; 674 return true; 675 } 676 #endif 677 678 /// \brief Finds the join node. 792 793 /// Finds the join node. 679 794 Node findJoinNode() { 680 // Finding the join node 681 Node u = graph.source(in_edge); 682 Node v = graph.target(in_edge); 795 Node u = _graph.source(_in_edge); 796 Node v = _graph.target(_in_edge); 683 797 while (u != v) { 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];798 if (_depth[u] == _depth[v]) { 799 u = _parent[u]; 800 v = _parent[v]; 801 } 802 else if (_depth[u] > _depth[v]) u = _parent[u]; 803 else v = _parent[v]; 690 804 } 691 805 return u; … … 697 811 // Initializing first and second nodes according to the direction 698 812 // of the cycle 699 if ( state[in_edge] ==LOWER) {700 first = graph.source(in_edge);701 second = graph.target(in_edge);813 if (_state[_in_edge] == STATE_LOWER) { 814 first = _graph.source(_in_edge); 815 second = _graph.target(_in_edge); 702 816 } else { 703 first = graph.target(in_edge);704 second = graph.source(in_edge);705 } 706 delta = capacity[in_edge];817 first = _graph.target(_in_edge); 818 second = _graph.source(_in_edge); 819 } 820 delta = _capacity[_in_edge]; 707 821 bool result = false; 708 822 Capacity d; … … 711 825 // Searching the cycle along the path form the first node to the 712 826 // root node 713 for (Node u = first; u != join; u = parent[u]) {714 e = pred_edge[u];715 d = forward[u] ? flow[e] : capacity[e] flow[e];827 for (Node u = first; u != join; u = _parent[u]) { 828 e = _pred_edge[u]; 829 d = _forward[u] ? _flow[e] : _capacity[e]  _flow[e]; 716 830 if (d < delta) { 717 831 delta = d; … … 724 838 // Searching the cycle along the path form the second node to the 725 839 // root node 726 for (Node u = second; u != join; u = parent[u]) {727 e = pred_edge[u];728 d = forward[u] ? capacity[e]  flow[e] :flow[e];840 for (Node u = second; u != join; u = _parent[u]) { 841 e = _pred_edge[u]; 842 d = _forward[u] ? _capacity[e]  _flow[e] : _flow[e]; 729 843 if (d <= delta) { 730 844 delta = d; … … 738 852 } 739 853 740 /// \briefChanges \c flow and \c state edge maps.854 /// Changes \c flow and \c state edge maps. 741 855 void changeFlows(bool change) { 742 856 // Augmenting along the cycle 743 857 if (delta > 0) { 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;858 Capacity val = _state[_in_edge] * delta; 859 _flow[_in_edge] += val; 860 for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) { 861 _flow[_pred_edge[u]] += _forward[u] ? val : val; 862 } 863 for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) { 864 _flow[_pred_edge[u]] += _forward[u] ? val : val; 751 865 } 752 866 } 753 867 // Updating the state of the entering and leaving edges 754 868 if (change) { 755 state[in_edge] =TREE;756 state[pred_edge[u_out]] =757 ( flow[pred_edge[u_out]] == 0) ? LOWER :UPPER;869 _state[_in_edge] = STATE_TREE; 870 _state[_pred_edge[u_out]] = 871 (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; 758 872 } else { 759 state[in_edge] = state[in_edge];760 } 761 } 762 763 /// \briefUpdates \c thread and \c parent node maps.873 _state[_in_edge] = _state[_in_edge]; 874 } 875 } 876 877 /// Updates \c thread and \c parent node maps. 764 878 void updateThreadParent() { 765 879 Node u; 766 v_out = parent[u_out];880 v_out = _parent[u_out]; 767 881 768 882 // Handling the case when join and v_out coincide 769 883 bool par_first = false; 770 884 if (join == v_out) { 771 for (u = join; u != u_in && u != v_in; u = thread[u]) ;885 for (u = join; u != u_in && u != v_in; u = _thread[u]) ; 772 886 if (u == v_in) { 773 887 par_first = true; 774 while ( thread[u] != u_out) u =thread[u];888 while (_thread[u] != u_out) u = _thread[u]; 775 889 first = u; 776 890 } … … 779 893 // Finding the last successor of u_in (u) and the node after it 780 894 // (right) according to the thread index 781 for (u = u_in; depth[thread[u]] > depth[u_in]; u =thread[u]) ;782 right = thread[u];783 if ( thread[v_in] == u_out) {784 for (last = u; depth[last] > depth[u_out]; last =thread[last]) ;785 if (last == u_out) last = thread[last];786 } 787 else last = thread[v_in];895 for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ; 896 right = _thread[u]; 897 if (_thread[v_in] == u_out) { 898 for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ; 899 if (last == u_out) last = _thread[last]; 900 } 901 else last = _thread[v_in]; 788 902 789 903 // Updating stem nodes 790 thread[v_in] = stem = u_in;904 _thread[v_in] = stem = u_in; 791 905 par_stem = v_in; 792 906 while (stem != u_out) { 793 thread[u] = new_stem =parent[stem];907 _thread[u] = new_stem = _parent[stem]; 794 908 795 909 // Finding the node just before the stem node (u) according to 796 910 // the original thread index 797 for (u = new_stem; thread[u] != stem; u =thread[u]) ;798 thread[u] = right;911 for (u = new_stem; _thread[u] != stem; u = _thread[u]) ; 912 _thread[u] = right; 799 913 800 914 // Changing the parent node of stem and shifting stem and 801 915 // par_stem nodes 802 parent[stem] = par_stem;916 _parent[stem] = par_stem; 803 917 par_stem = stem; 804 918 stem = new_stem; … … 806 920 // Finding the last successor of stem (u) and the node after it 807 921 // (right) according to the thread index 808 for (u = stem; depth[thread[u]] > depth[stem]; u =thread[u]) ;809 right = thread[u];810 } 811 parent[u_out] = par_stem;812 thread[u] = last;922 for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ; 923 right = _thread[u]; 924 } 925 _parent[u_out] = par_stem; 926 _thread[u] = last; 813 927 814 928 if (join == v_out && par_first) { 815 if (first != v_in) thread[first] = right;929 if (first != v_in) _thread[first] = right; 816 930 } else { 817 for (u = v_out; thread[u] != u_out; u =thread[u]) ;818 thread[u] = right;819 } 820 } 821 822 /// \briefUpdates \c pred_edge and \c forward node maps.931 for (u = v_out; _thread[u] != u_out; u = _thread[u]) ; 932 _thread[u] = right; 933 } 934 } 935 936 /// Updates \c pred_edge and \c forward node maps. 823 937 void updatePredEdge() { 824 938 Node u = u_out, v; 825 939 while (u != u_in) { 826 v = parent[u];827 pred_edge[u] =pred_edge[v];828 forward[u] = !forward[v];940 v = _parent[u]; 941 _pred_edge[u] = _pred_edge[v]; 942 _forward[u] = !_forward[v]; 829 943 u = v; 830 944 } 831 pred_edge[u_in] =in_edge;832 forward[u_in] = (u_in == graph.source(in_edge));833 } 834 835 /// \briefUpdates \c depth and \c potential node maps.945 _pred_edge[u_in] = _in_edge; 946 _forward[u_in] = (u_in == _graph.source(_in_edge)); 947 } 948 949 /// Updates \c depth and \c potential node maps. 836 950 void updateDepthPotential() { 837 depth[u_in] =depth[v_in] + 1;838 potential[u_in] =forward[u_in] ?839 potential[v_in] + cost[pred_edge[u_in]] :840 potential[v_in]  cost[pred_edge[u_in]];841 842 Node u = thread[u_in], v;951 _depth[u_in] = _depth[v_in] + 1; 952 _potential[u_in] = _forward[u_in] ? 953 _potential[v_in]  _cost[_pred_edge[u_in]] : 954 _potential[v_in] + _cost[_pred_edge[u_in]]; 955 956 Node u = _thread[u_in], v; 843 957 while (true) { 844 v = parent[u];958 v = _parent[u]; 845 959 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]; 852 } 853 } 854 855 /// \brief Executes the algorithm. 960 _depth[u] = _depth[v] + 1; 961 _potential[u] = _forward[u] ? 962 _potential[v]  _cost[_pred_edge[u]] : 963 _potential[v] + _cost[_pred_edge[u]]; 964 if (_depth[u] <= _depth[v_in]) break; 965 u = _thread[u]; 966 } 967 } 968 969 /// Executes the algorithm. 970 bool start(PivotRuleEnum pivot_rule) { 971 switch (pivot_rule) { 972 case FIRST_ELIGIBLE_PIVOT: 973 return start<FirstEligiblePivotRule>(); 974 case BEST_ELIGIBLE_PIVOT: 975 return start<BestEligiblePivotRule>(); 976 case BLOCK_SEARCH_PIVOT: 977 return start<BlockSearchPivotRule>(); 978 case LIMITED_SEARCH_PIVOT: 979 return start<LimitedSearchPivotRule>(); 980 case CANDIDATE_LIST_PIVOT: 981 return start<CandidateListPivotRule>(); 982 case COMBINED_PIVOT: 983 if ( countEdges(_graph) / countNodes(_graph) <= 984 COMBINED_PIVOT_MAX_DEG ) 985 return start<LimitedSearchPivotRule>(); 986 else 987 return start<BlockSearchPivotRule>(); 988 } 989 return false; 990 } 991 992 template<class PivotRuleImplementation> 856 993 bool start() { 857 // Processing pivots 858 #ifdef _DEBUG_ITER_ 859 int iter_num = 0; 860 #endif 861 #if defined(FIRST_ELIGIBLE_PIVOT)  defined(EDGE_BLOCK_PIVOT) 862 EdgeIt next_edge(graph); 863 while (findEnteringEdge(next_edge)) 864 #else 865 while (findEnteringEdge()) 866 #endif 867 { 994 PivotRuleImplementation pivot(*this); 995 996 // Executing the network simplex algorithm 997 while (pivot.findEnteringEdge()) { 868 998 join = findJoinNode(); 869 999 bool change = findLeavingEdge(); … … 874 1004 updateDepthPotential(); 875 1005 } 876 #ifdef _DEBUG_ITER_ 877 ++iter_num; 878 #endif 879 } 880 881 #ifdef _DEBUG_ITER_ 882 std::cout << "Network Simplex algorithm finished. " << iter_num 883 << " pivot iterations performed." << std::endl; 884 #endif 885 886 // Checking if the flow amount equals zero on all the 887 // artificial edges 888 for (InEdgeIt e(graph, root); e != INVALID; ++e) 889 if (flow[e] > 0) return false; 890 for (OutEdgeIt e(graph, root); e != INVALID; ++e) 891 if (flow[e] > 0) return false; 892 893 // Copying flow values to flow_result 894 if (lower) { 895 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) 896 flow_result[e] = (*lower)[e] + flow[edge_ref[e]]; 1006 } 1007 1008 // Checking if the flow amount equals zero on all the artificial 1009 // edges 1010 for (InEdgeIt e(_graph, _root); e != INVALID; ++e) 1011 if (_flow[e] > 0) return false; 1012 for (OutEdgeIt e(_graph, _root); e != INVALID; ++e) 1013 if (_flow[e] > 0) return false; 1014 1015 // Copying flow values to _flow_result 1016 if (_lower) { 1017 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) 1018 _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]]; 897 1019 } else { 898 for (typename Graph::EdgeIt e( graph_ref); e != INVALID; ++e)899 flow_result[e] = flow[edge_ref[e]];900 } 901 // Copying potential values to potential_result902 for (typename Graph::NodeIt n( graph_ref); n != INVALID; ++n)903 potential_result[n] = potential[node_ref[n]];1020 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) 1021 _flow_result[e] = _flow[_edge_ref[e]]; 1022 } 1023 // Copying potential values to _potential_result 1024 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) 1025 _potential_result[n] = _potential[_node_ref[n]]; 904 1026 905 1027 return true;
Note: See TracChangeset
for help on using the changeset viewer.