[2440] | 1 | /* -*- C++ -*- |
---|
| 2 | * |
---|
| 3 | * This file is a part of LEMON, a generic C++ optimization library |
---|
| 4 | * |
---|
[2553] | 5 | * Copyright (C) 2003-2008 |
---|
[2440] | 6 | * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
---|
| 7 | * (Egervary Research Group on Combinatorial Optimization, EGRES). |
---|
| 8 | * |
---|
| 9 | * Permission to use, modify and distribute this software is granted |
---|
| 10 | * provided that this copyright notice appears in all copies. For |
---|
| 11 | * precise terms see the accompanying LICENSE file. |
---|
| 12 | * |
---|
| 13 | * This software is provided "AS IS" with no warranty of any kind, |
---|
| 14 | * express or implied, and with no claim as to its suitability for any |
---|
| 15 | * purpose. |
---|
| 16 | * |
---|
| 17 | */ |
---|
| 18 | |
---|
| 19 | #ifndef LEMON_NETWORK_SIMPLEX_H |
---|
| 20 | #define LEMON_NETWORK_SIMPLEX_H |
---|
| 21 | |
---|
| 22 | /// \ingroup min_cost_flow |
---|
| 23 | /// |
---|
| 24 | /// \file |
---|
[2575] | 25 | /// \brief Network simplex algorithm for finding a minimum cost flow. |
---|
[2440] | 26 | |
---|
[2575] | 27 | #include <vector> |
---|
[2440] | 28 | #include <limits> |
---|
[2575] | 29 | |
---|
[2509] | 30 | #include <lemon/graph_adaptor.h> |
---|
| 31 | #include <lemon/graph_utils.h> |
---|
[2440] | 32 | #include <lemon/smart_graph.h> |
---|
[2575] | 33 | #include <lemon/math.h> |
---|
[2440] | 34 | |
---|
| 35 | namespace lemon { |
---|
| 36 | |
---|
| 37 | /// \addtogroup min_cost_flow |
---|
| 38 | /// @{ |
---|
| 39 | |
---|
[2619] | 40 | /// \brief Implementation of the primal network simplex algorithm |
---|
| 41 | /// for finding a minimum cost flow. |
---|
[2440] | 42 | /// |
---|
[2619] | 43 | /// \ref NetworkSimplex implements the primal network simplex algorithm |
---|
| 44 | /// for finding a minimum cost flow. |
---|
[2440] | 45 | /// |
---|
[2575] | 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. |
---|
[2440] | 51 | /// |
---|
| 52 | /// \warning |
---|
[2575] | 53 | /// - Edge capacities and costs should be \e non-negative \e integers. |
---|
| 54 | /// - Supply values should be \e signed \e integers. |
---|
[2581] | 55 | /// - The value types of the maps should be convertible to each other. |
---|
| 56 | /// - \c CostMap::Value must be signed type. |
---|
[2575] | 57 | /// |
---|
[2619] | 58 | /// \note \ref NetworkSimplex provides five different pivot rule |
---|
[2575] | 59 | /// implementations that significantly affect the efficiency of the |
---|
| 60 | /// algorithm. |
---|
[2619] | 61 | /// By default "Block Search" pivot rule is used, which proved to be |
---|
| 62 | /// by far the most efficient according to our benchmark tests. |
---|
| 63 | /// However another pivot rule can be selected using \ref run() |
---|
| 64 | /// function with the proper parameter. |
---|
[2440] | 65 | /// |
---|
| 66 | /// \author Peter Kovacs |
---|
[2533] | 67 | template < typename Graph, |
---|
| 68 | typename LowerMap = typename Graph::template EdgeMap<int>, |
---|
[2575] | 69 | typename CapacityMap = typename Graph::template EdgeMap<int>, |
---|
[2533] | 70 | typename CostMap = typename Graph::template EdgeMap<int>, |
---|
[2575] | 71 | typename SupplyMap = typename Graph::template NodeMap<int> > |
---|
[2440] | 72 | class NetworkSimplex |
---|
| 73 | { |
---|
| 74 | typedef typename CapacityMap::Value Capacity; |
---|
| 75 | typedef typename CostMap::Value Cost; |
---|
| 76 | typedef typename SupplyMap::Value Supply; |
---|
| 77 | |
---|
| 78 | typedef SmartGraph SGraph; |
---|
[2556] | 79 | GRAPH_TYPEDEFS(typename SGraph); |
---|
[2440] | 80 | |
---|
| 81 | typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap; |
---|
| 82 | typedef typename SGraph::template EdgeMap<Cost> SCostMap; |
---|
| 83 | typedef typename SGraph::template NodeMap<Supply> SSupplyMap; |
---|
| 84 | typedef typename SGraph::template NodeMap<Cost> SPotentialMap; |
---|
| 85 | |
---|
| 86 | typedef typename SGraph::template NodeMap<int> IntNodeMap; |
---|
| 87 | typedef typename SGraph::template NodeMap<bool> BoolNodeMap; |
---|
| 88 | typedef typename SGraph::template NodeMap<Node> NodeNodeMap; |
---|
| 89 | typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap; |
---|
| 90 | typedef typename SGraph::template EdgeMap<int> IntEdgeMap; |
---|
[2619] | 91 | typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap; |
---|
[2440] | 92 | |
---|
| 93 | typedef typename Graph::template NodeMap<Node> NodeRefMap; |
---|
| 94 | typedef typename Graph::template EdgeMap<Edge> EdgeRefMap; |
---|
| 95 | |
---|
[2619] | 96 | typedef std::vector<Edge> EdgeVector; |
---|
| 97 | |
---|
[2440] | 98 | public: |
---|
| 99 | |
---|
[2556] | 100 | /// The type of the flow map. |
---|
[2440] | 101 | typedef typename Graph::template EdgeMap<Capacity> FlowMap; |
---|
[2556] | 102 | /// The type of the potential map. |
---|
[2440] | 103 | typedef typename Graph::template NodeMap<Cost> PotentialMap; |
---|
| 104 | |
---|
[2575] | 105 | public: |
---|
[2440] | 106 | |
---|
[2575] | 107 | /// Enum type to select the pivot rule used by \ref run(). |
---|
| 108 | enum PivotRuleEnum { |
---|
| 109 | FIRST_ELIGIBLE_PIVOT, |
---|
| 110 | BEST_ELIGIBLE_PIVOT, |
---|
| 111 | BLOCK_SEARCH_PIVOT, |
---|
| 112 | CANDIDATE_LIST_PIVOT, |
---|
[2619] | 113 | ALTERING_LIST_PIVOT |
---|
[2575] | 114 | }; |
---|
| 115 | |
---|
| 116 | private: |
---|
| 117 | |
---|
| 118 | /// \brief Map adaptor class for handling reduced edge costs. |
---|
| 119 | /// |
---|
[2556] | 120 | /// Map adaptor class for handling reduced edge costs. |
---|
[2440] | 121 | class ReducedCostMap : public MapBase<Edge, Cost> |
---|
| 122 | { |
---|
| 123 | private: |
---|
| 124 | |
---|
[2575] | 125 | const SGraph &_gr; |
---|
| 126 | const SCostMap &_cost_map; |
---|
| 127 | const SPotentialMap &_pot_map; |
---|
[2440] | 128 | |
---|
| 129 | public: |
---|
| 130 | |
---|
[2575] | 131 | ///\e |
---|
| 132 | ReducedCostMap( const SGraph &gr, |
---|
| 133 | const SCostMap &cost_map, |
---|
| 134 | const SPotentialMap &pot_map ) : |
---|
[2579] | 135 | _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {} |
---|
[2440] | 136 | |
---|
[2575] | 137 | ///\e |
---|
[2509] | 138 | Cost operator[](const Edge &e) const { |
---|
[2575] | 139 | return _cost_map[e] + _pot_map[_gr.source(e)] |
---|
| 140 | - _pot_map[_gr.target(e)]; |
---|
[2440] | 141 | } |
---|
| 142 | |
---|
| 143 | }; //class ReducedCostMap |
---|
| 144 | |
---|
[2575] | 145 | private: |
---|
[2440] | 146 | |
---|
[2575] | 147 | /// \brief Implementation of the "First Eligible" pivot rule for the |
---|
| 148 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 149 | /// |
---|
| 150 | /// This class implements the "First Eligible" pivot rule |
---|
| 151 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
[2619] | 152 | /// |
---|
| 153 | /// For more information see \ref NetworkSimplex::run(). |
---|
[2575] | 154 | class FirstEligiblePivotRule |
---|
| 155 | { |
---|
| 156 | private: |
---|
[2440] | 157 | |
---|
[2619] | 158 | // References to the NetworkSimplex class |
---|
[2575] | 159 | NetworkSimplex &_ns; |
---|
[2619] | 160 | EdgeVector &_edges; |
---|
| 161 | |
---|
| 162 | int _next_edge; |
---|
[2440] | 163 | |
---|
[2575] | 164 | public: |
---|
[2440] | 165 | |
---|
[2619] | 166 | /// Constructor |
---|
| 167 | FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) : |
---|
| 168 | _ns(ns), _edges(edges), _next_edge(0) {} |
---|
[2575] | 169 | |
---|
[2619] | 170 | /// Find next entering edge |
---|
| 171 | inline bool findEnteringEdge() { |
---|
| 172 | Edge e; |
---|
| 173 | for (int i = _next_edge; i < int(_edges.size()); ++i) { |
---|
| 174 | e = _edges[i]; |
---|
[2575] | 175 | if (_ns._state[e] * _ns._red_cost[e] < 0) { |
---|
| 176 | _ns._in_edge = e; |
---|
[2619] | 177 | _next_edge = i + 1; |
---|
[2575] | 178 | return true; |
---|
| 179 | } |
---|
| 180 | } |
---|
[2619] | 181 | for (int i = 0; i < _next_edge; ++i) { |
---|
| 182 | e = _edges[i]; |
---|
[2575] | 183 | if (_ns._state[e] * _ns._red_cost[e] < 0) { |
---|
| 184 | _ns._in_edge = e; |
---|
[2619] | 185 | _next_edge = i + 1; |
---|
[2575] | 186 | return true; |
---|
| 187 | } |
---|
| 188 | } |
---|
| 189 | return false; |
---|
| 190 | } |
---|
| 191 | }; //class FirstEligiblePivotRule |
---|
| 192 | |
---|
| 193 | /// \brief Implementation of the "Best Eligible" pivot rule for the |
---|
| 194 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 195 | /// |
---|
| 196 | /// This class implements the "Best Eligible" pivot rule |
---|
| 197 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
[2619] | 198 | /// |
---|
| 199 | /// For more information see \ref NetworkSimplex::run(). |
---|
[2575] | 200 | class BestEligiblePivotRule |
---|
| 201 | { |
---|
| 202 | private: |
---|
| 203 | |
---|
[2619] | 204 | // References to the NetworkSimplex class |
---|
[2575] | 205 | NetworkSimplex &_ns; |
---|
[2619] | 206 | EdgeVector &_edges; |
---|
[2575] | 207 | |
---|
| 208 | public: |
---|
| 209 | |
---|
[2619] | 210 | /// Constructor |
---|
| 211 | BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) : |
---|
| 212 | _ns(ns), _edges(edges) {} |
---|
[2575] | 213 | |
---|
[2619] | 214 | /// Find next entering edge |
---|
| 215 | inline bool findEnteringEdge() { |
---|
[2575] | 216 | Cost min = 0; |
---|
[2619] | 217 | Edge e; |
---|
| 218 | for (int i = 0; i < int(_edges.size()); ++i) { |
---|
| 219 | e = _edges[i]; |
---|
[2575] | 220 | if (_ns._state[e] * _ns._red_cost[e] < min) { |
---|
| 221 | min = _ns._state[e] * _ns._red_cost[e]; |
---|
| 222 | _ns._in_edge = e; |
---|
| 223 | } |
---|
| 224 | } |
---|
| 225 | return min < 0; |
---|
| 226 | } |
---|
| 227 | }; //class BestEligiblePivotRule |
---|
| 228 | |
---|
| 229 | /// \brief Implementation of the "Block Search" pivot rule for the |
---|
| 230 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 231 | /// |
---|
| 232 | /// This class implements the "Block Search" pivot rule |
---|
| 233 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
[2619] | 234 | /// |
---|
| 235 | /// For more information see \ref NetworkSimplex::run(). |
---|
[2575] | 236 | class BlockSearchPivotRule |
---|
| 237 | { |
---|
| 238 | private: |
---|
| 239 | |
---|
[2619] | 240 | // References to the NetworkSimplex class |
---|
[2575] | 241 | NetworkSimplex &_ns; |
---|
[2619] | 242 | EdgeVector &_edges; |
---|
| 243 | |
---|
[2575] | 244 | int _block_size; |
---|
[2619] | 245 | int _next_edge, _min_edge; |
---|
[2575] | 246 | |
---|
| 247 | public: |
---|
| 248 | |
---|
[2619] | 249 | /// Constructor |
---|
| 250 | BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) : |
---|
| 251 | _ns(ns), _edges(edges), _next_edge(0), _min_edge(0) |
---|
[2575] | 252 | { |
---|
[2619] | 253 | // The main parameters of the pivot rule |
---|
| 254 | const double BLOCK_SIZE_FACTOR = 2.0; |
---|
| 255 | const int MIN_BLOCK_SIZE = 10; |
---|
| 256 | |
---|
| 257 | _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())), |
---|
| 258 | MIN_BLOCK_SIZE ); |
---|
[2575] | 259 | } |
---|
| 260 | |
---|
[2619] | 261 | /// Find next entering edge |
---|
| 262 | inline bool findEnteringEdge() { |
---|
[2575] | 263 | Cost curr, min = 0; |
---|
[2619] | 264 | Edge e; |
---|
| 265 | int cnt = _block_size; |
---|
| 266 | int i; |
---|
| 267 | for (i = _next_edge; i < int(_edges.size()); ++i) { |
---|
| 268 | e = _edges[i]; |
---|
[2575] | 269 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
---|
| 270 | min = curr; |
---|
[2619] | 271 | _min_edge = i; |
---|
[2575] | 272 | } |
---|
[2619] | 273 | if (--cnt == 0) { |
---|
[2575] | 274 | if (min < 0) break; |
---|
[2619] | 275 | cnt = _block_size; |
---|
[2575] | 276 | } |
---|
| 277 | } |
---|
[2619] | 278 | if (min == 0 || cnt > 0) { |
---|
| 279 | for (i = 0; i < _next_edge; ++i) { |
---|
| 280 | e = _edges[i]; |
---|
[2575] | 281 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
---|
| 282 | min = curr; |
---|
[2619] | 283 | _min_edge = i; |
---|
[2575] | 284 | } |
---|
[2619] | 285 | if (--cnt == 0) { |
---|
[2575] | 286 | if (min < 0) break; |
---|
[2619] | 287 | cnt = _block_size; |
---|
[2575] | 288 | } |
---|
| 289 | } |
---|
| 290 | } |
---|
[2619] | 291 | if (min >= 0) return false; |
---|
| 292 | _ns._in_edge = _edges[_min_edge]; |
---|
| 293 | _next_edge = i; |
---|
| 294 | return true; |
---|
[2575] | 295 | } |
---|
| 296 | }; //class BlockSearchPivotRule |
---|
| 297 | |
---|
| 298 | /// \brief Implementation of the "Candidate List" pivot rule for the |
---|
| 299 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 300 | /// |
---|
| 301 | /// This class implements the "Candidate List" pivot rule |
---|
| 302 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
[2619] | 303 | /// |
---|
| 304 | /// For more information see \ref NetworkSimplex::run(). |
---|
[2575] | 305 | class CandidateListPivotRule |
---|
| 306 | { |
---|
| 307 | private: |
---|
| 308 | |
---|
[2619] | 309 | // References to the NetworkSimplex class |
---|
[2575] | 310 | NetworkSimplex &_ns; |
---|
[2619] | 311 | EdgeVector &_edges; |
---|
[2575] | 312 | |
---|
[2619] | 313 | EdgeVector _candidates; |
---|
| 314 | int _list_length, _minor_limit; |
---|
| 315 | int _curr_length, _minor_count; |
---|
| 316 | int _next_edge, _min_edge; |
---|
[2575] | 317 | |
---|
| 318 | public: |
---|
| 319 | |
---|
[2619] | 320 | /// Constructor |
---|
| 321 | CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) : |
---|
| 322 | _ns(ns), _edges(edges), _next_edge(0), _min_edge(0) |
---|
[2575] | 323 | { |
---|
[2619] | 324 | // The main parameters of the pivot rule |
---|
| 325 | const double LIST_LENGTH_FACTOR = 1.0; |
---|
| 326 | const int MIN_LIST_LENGTH = 10; |
---|
| 327 | const double MINOR_LIMIT_FACTOR = 0.1; |
---|
| 328 | const int MIN_MINOR_LIMIT = 3; |
---|
| 329 | |
---|
| 330 | _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())), |
---|
| 331 | MIN_LIST_LENGTH ); |
---|
| 332 | _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), |
---|
| 333 | MIN_MINOR_LIMIT ); |
---|
| 334 | _curr_length = _minor_count = 0; |
---|
| 335 | _candidates.resize(_list_length); |
---|
[2575] | 336 | } |
---|
| 337 | |
---|
[2619] | 338 | /// Find next entering edge |
---|
| 339 | inline bool findEnteringEdge() { |
---|
[2575] | 340 | Cost min, curr; |
---|
[2619] | 341 | if (_curr_length > 0 && _minor_count < _minor_limit) { |
---|
| 342 | // Minor iteration: selecting the best eligible edge from |
---|
| 343 | // the current candidate list |
---|
[2575] | 344 | ++_minor_count; |
---|
| 345 | Edge e; |
---|
| 346 | min = 0; |
---|
[2619] | 347 | for (int i = 0; i < _curr_length; ++i) { |
---|
[2575] | 348 | e = _candidates[i]; |
---|
[2619] | 349 | curr = _ns._state[e] * _ns._red_cost[e]; |
---|
| 350 | if (curr < min) { |
---|
[2575] | 351 | min = curr; |
---|
| 352 | _ns._in_edge = e; |
---|
| 353 | } |
---|
[2619] | 354 | if (curr >= 0) { |
---|
| 355 | _candidates[i--] = _candidates[--_curr_length]; |
---|
| 356 | } |
---|
[2575] | 357 | } |
---|
| 358 | if (min < 0) return true; |
---|
| 359 | } |
---|
| 360 | |
---|
[2619] | 361 | // Major iteration: building a new candidate list |
---|
| 362 | Edge e; |
---|
[2575] | 363 | min = 0; |
---|
[2619] | 364 | _curr_length = 0; |
---|
| 365 | int i; |
---|
| 366 | for (i = _next_edge; i < int(_edges.size()); ++i) { |
---|
| 367 | e = _edges[i]; |
---|
[2575] | 368 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { |
---|
[2619] | 369 | _candidates[_curr_length++] = e; |
---|
[2575] | 370 | if (curr < min) { |
---|
| 371 | min = curr; |
---|
[2619] | 372 | _min_edge = i; |
---|
[2575] | 373 | } |
---|
[2619] | 374 | if (_curr_length == _list_length) break; |
---|
[2575] | 375 | } |
---|
| 376 | } |
---|
[2619] | 377 | if (_curr_length < _list_length) { |
---|
| 378 | for (i = 0; i < _next_edge; ++i) { |
---|
| 379 | e = _edges[i]; |
---|
[2575] | 380 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { |
---|
[2619] | 381 | _candidates[_curr_length++] = e; |
---|
[2575] | 382 | if (curr < min) { |
---|
| 383 | min = curr; |
---|
[2619] | 384 | _min_edge = i; |
---|
[2575] | 385 | } |
---|
[2619] | 386 | if (_curr_length == _list_length) break; |
---|
[2575] | 387 | } |
---|
| 388 | } |
---|
| 389 | } |
---|
[2619] | 390 | if (_curr_length == 0) return false; |
---|
[2575] | 391 | _minor_count = 1; |
---|
[2619] | 392 | _ns._in_edge = _edges[_min_edge]; |
---|
| 393 | _next_edge = i; |
---|
[2575] | 394 | return true; |
---|
| 395 | } |
---|
| 396 | }; //class CandidateListPivotRule |
---|
| 397 | |
---|
[2619] | 398 | /// \brief Implementation of the "Altering Candidate List" pivot rule |
---|
| 399 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
| 400 | /// |
---|
| 401 | /// This class implements the "Altering Candidate List" pivot rule |
---|
| 402 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
| 403 | /// |
---|
| 404 | /// For more information see \ref NetworkSimplex::run(). |
---|
| 405 | class AlteringListPivotRule |
---|
| 406 | { |
---|
| 407 | private: |
---|
| 408 | |
---|
| 409 | // References to the NetworkSimplex class |
---|
| 410 | NetworkSimplex &_ns; |
---|
| 411 | EdgeVector &_edges; |
---|
| 412 | |
---|
| 413 | EdgeVector _candidates; |
---|
| 414 | SCostMap _cand_cost; |
---|
| 415 | int _block_size, _head_length, _curr_length; |
---|
| 416 | int _next_edge; |
---|
| 417 | |
---|
| 418 | // Functor class to compare edges during sort of the candidate list |
---|
| 419 | class SortFunc |
---|
| 420 | { |
---|
| 421 | private: |
---|
| 422 | const SCostMap &_map; |
---|
| 423 | public: |
---|
| 424 | SortFunc(const SCostMap &map) : _map(map) {} |
---|
| 425 | bool operator()(const Edge &e1, const Edge &e2) { |
---|
| 426 | return _map[e1] < _map[e2]; |
---|
| 427 | } |
---|
| 428 | }; |
---|
| 429 | |
---|
| 430 | SortFunc _sort_func; |
---|
| 431 | |
---|
| 432 | public: |
---|
| 433 | |
---|
| 434 | /// Constructor |
---|
| 435 | AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) : |
---|
| 436 | _ns(ns), _edges(edges), _cand_cost(_ns._graph), |
---|
| 437 | _next_edge(0), _sort_func(_cand_cost) |
---|
| 438 | { |
---|
| 439 | // The main parameters of the pivot rule |
---|
| 440 | const double BLOCK_SIZE_FACTOR = 1.0; |
---|
| 441 | const int MIN_BLOCK_SIZE = 10; |
---|
| 442 | const double HEAD_LENGTH_FACTOR = 0.1; |
---|
| 443 | const int MIN_HEAD_LENGTH = 5; |
---|
| 444 | |
---|
| 445 | _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())), |
---|
| 446 | MIN_BLOCK_SIZE ); |
---|
| 447 | _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), |
---|
| 448 | MIN_HEAD_LENGTH ); |
---|
| 449 | _candidates.resize(_head_length + _block_size); |
---|
| 450 | _curr_length = 0; |
---|
| 451 | } |
---|
| 452 | |
---|
| 453 | /// Find next entering edge |
---|
| 454 | inline bool findEnteringEdge() { |
---|
| 455 | // Checking the current candidate list |
---|
| 456 | Edge e; |
---|
| 457 | for (int idx = 0; idx < _curr_length; ++idx) { |
---|
| 458 | e = _candidates[idx]; |
---|
| 459 | if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) { |
---|
| 460 | _candidates[idx--] = _candidates[--_curr_length]; |
---|
| 461 | } |
---|
| 462 | } |
---|
| 463 | |
---|
| 464 | // Extending the list |
---|
| 465 | int cnt = _block_size; |
---|
| 466 | int last_edge = 0; |
---|
| 467 | int limit = _head_length; |
---|
| 468 | for (int i = _next_edge; i < int(_edges.size()); ++i) { |
---|
| 469 | e = _edges[i]; |
---|
| 470 | if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) { |
---|
| 471 | _candidates[_curr_length++] = e; |
---|
| 472 | last_edge = i; |
---|
| 473 | } |
---|
| 474 | if (--cnt == 0) { |
---|
| 475 | if (_curr_length > limit) break; |
---|
| 476 | limit = 0; |
---|
| 477 | cnt = _block_size; |
---|
| 478 | } |
---|
| 479 | } |
---|
| 480 | if (_curr_length <= limit) { |
---|
| 481 | for (int i = 0; i < _next_edge; ++i) { |
---|
| 482 | e = _edges[i]; |
---|
| 483 | if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) { |
---|
| 484 | _candidates[_curr_length++] = e; |
---|
| 485 | last_edge = i; |
---|
| 486 | } |
---|
| 487 | if (--cnt == 0) { |
---|
| 488 | if (_curr_length > limit) break; |
---|
| 489 | limit = 0; |
---|
| 490 | cnt = _block_size; |
---|
| 491 | } |
---|
| 492 | } |
---|
| 493 | } |
---|
| 494 | if (_curr_length == 0) return false; |
---|
| 495 | _next_edge = last_edge + 1; |
---|
| 496 | |
---|
| 497 | // Sorting the list partially |
---|
| 498 | EdgeVector::iterator sort_end = _candidates.begin(); |
---|
| 499 | EdgeVector::iterator vector_end = _candidates.begin(); |
---|
| 500 | for (int idx = 0; idx < _curr_length; ++idx) { |
---|
| 501 | ++vector_end; |
---|
| 502 | if (idx <= _head_length) ++sort_end; |
---|
| 503 | } |
---|
| 504 | partial_sort(_candidates.begin(), sort_end, vector_end, _sort_func); |
---|
| 505 | |
---|
| 506 | _ns._in_edge = _candidates[0]; |
---|
| 507 | if (_curr_length > _head_length) { |
---|
| 508 | _candidates[0] = _candidates[_head_length - 1]; |
---|
| 509 | _curr_length = _head_length - 1; |
---|
| 510 | } else { |
---|
| 511 | _candidates[0] = _candidates[_curr_length - 1]; |
---|
| 512 | --_curr_length; |
---|
| 513 | } |
---|
| 514 | |
---|
| 515 | return true; |
---|
| 516 | } |
---|
| 517 | }; //class AlteringListPivotRule |
---|
| 518 | |
---|
[2575] | 519 | private: |
---|
| 520 | |
---|
[2579] | 521 | // State constants for edges |
---|
| 522 | enum EdgeStateEnum { |
---|
| 523 | STATE_UPPER = -1, |
---|
| 524 | STATE_TREE = 0, |
---|
| 525 | STATE_LOWER = 1 |
---|
| 526 | }; |
---|
[2575] | 527 | |
---|
| 528 | private: |
---|
| 529 | |
---|
| 530 | // The directed graph the algorithm runs on |
---|
| 531 | SGraph _graph; |
---|
| 532 | // The original graph |
---|
| 533 | const Graph &_graph_ref; |
---|
| 534 | // The original lower bound map |
---|
| 535 | const LowerMap *_lower; |
---|
| 536 | // The capacity map |
---|
| 537 | SCapacityMap _capacity; |
---|
| 538 | // The cost map |
---|
| 539 | SCostMap _cost; |
---|
| 540 | // The supply map |
---|
| 541 | SSupplyMap _supply; |
---|
| 542 | bool _valid_supply; |
---|
| 543 | |
---|
| 544 | // Edge map of the current flow |
---|
| 545 | SCapacityMap _flow; |
---|
| 546 | // Node map of the current potentials |
---|
| 547 | SPotentialMap _potential; |
---|
| 548 | |
---|
| 549 | // The depth node map of the spanning tree structure |
---|
| 550 | IntNodeMap _depth; |
---|
| 551 | // The parent node map of the spanning tree structure |
---|
| 552 | NodeNodeMap _parent; |
---|
| 553 | // The pred_edge node map of the spanning tree structure |
---|
| 554 | EdgeNodeMap _pred_edge; |
---|
| 555 | // The thread node map of the spanning tree structure |
---|
| 556 | NodeNodeMap _thread; |
---|
| 557 | // The forward node map of the spanning tree structure |
---|
| 558 | BoolNodeMap _forward; |
---|
| 559 | // The state edge map |
---|
| 560 | IntEdgeMap _state; |
---|
| 561 | // The root node of the starting spanning tree |
---|
| 562 | Node _root; |
---|
| 563 | |
---|
| 564 | // The reduced cost map |
---|
| 565 | ReducedCostMap _red_cost; |
---|
| 566 | |
---|
[2619] | 567 | // The non-artifical edges |
---|
| 568 | EdgeVector _edges; |
---|
| 569 | |
---|
[2575] | 570 | // Members for handling the original graph |
---|
[2581] | 571 | FlowMap *_flow_result; |
---|
| 572 | PotentialMap *_potential_result; |
---|
| 573 | bool _local_flow; |
---|
| 574 | bool _local_potential; |
---|
[2575] | 575 | NodeRefMap _node_ref; |
---|
| 576 | EdgeRefMap _edge_ref; |
---|
[2440] | 577 | |
---|
[2556] | 578 | // The entering edge of the current pivot iteration. |
---|
[2575] | 579 | Edge _in_edge; |
---|
| 580 | |
---|
[2556] | 581 | // Temporary nodes used in the current pivot iteration. |
---|
| 582 | Node join, u_in, v_in, u_out, v_out; |
---|
| 583 | Node right, first, second, last; |
---|
| 584 | Node stem, par_stem, new_stem; |
---|
| 585 | // The maximum augment amount along the found cycle in the current |
---|
| 586 | // pivot iteration. |
---|
| 587 | Capacity delta; |
---|
[2440] | 588 | |
---|
| 589 | public : |
---|
| 590 | |
---|
[2581] | 591 | /// \brief General constructor (with lower bounds). |
---|
[2440] | 592 | /// |
---|
[2581] | 593 | /// General constructor (with lower bounds). |
---|
[2440] | 594 | /// |
---|
[2575] | 595 | /// \param graph The directed graph the algorithm runs on. |
---|
| 596 | /// \param lower The lower bounds of the edges. |
---|
| 597 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 598 | /// \param cost The cost (length) values of the edges. |
---|
| 599 | /// \param supply The supply values of the nodes (signed). |
---|
| 600 | NetworkSimplex( const Graph &graph, |
---|
| 601 | const LowerMap &lower, |
---|
| 602 | const CapacityMap &capacity, |
---|
| 603 | const CostMap &cost, |
---|
| 604 | const SupplyMap &supply ) : |
---|
| 605 | _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), |
---|
| 606 | _cost(_graph), _supply(_graph), _flow(_graph), |
---|
| 607 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 608 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 609 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
[2623] | 610 | _flow_result(NULL), _potential_result(NULL), |
---|
[2581] | 611 | _local_flow(false), _local_potential(false), |
---|
[2575] | 612 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 613 | { |
---|
| 614 | // Checking the sum of supply values |
---|
| 615 | Supply sum = 0; |
---|
[2575] | 616 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
---|
| 617 | sum += supply[n]; |
---|
| 618 | if (!(_valid_supply = sum == 0)) return; |
---|
[2440] | 619 | |
---|
[2575] | 620 | // Copying _graph_ref to _graph |
---|
| 621 | _graph.reserveNode(countNodes(_graph_ref) + 1); |
---|
| 622 | _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); |
---|
| 623 | copyGraph(_graph, _graph_ref) |
---|
| 624 | .edgeMap(_cost, cost) |
---|
| 625 | .nodeRef(_node_ref) |
---|
| 626 | .edgeRef(_edge_ref) |
---|
[2556] | 627 | .run(); |
---|
[2440] | 628 | |
---|
[2556] | 629 | // Removing non-zero lower bounds |
---|
[2575] | 630 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { |
---|
| 631 | _capacity[_edge_ref[e]] = capacity[e] - lower[e]; |
---|
[2440] | 632 | } |
---|
[2575] | 633 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { |
---|
| 634 | Supply s = supply[n]; |
---|
| 635 | for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 636 | s += lower[e]; |
---|
| 637 | for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 638 | s -= lower[e]; |
---|
| 639 | _supply[_node_ref[n]] = s; |
---|
[2440] | 640 | } |
---|
| 641 | } |
---|
| 642 | |
---|
[2581] | 643 | /// \brief General constructor (without lower bounds). |
---|
[2440] | 644 | /// |
---|
[2581] | 645 | /// General constructor (without lower bounds). |
---|
[2440] | 646 | /// |
---|
[2575] | 647 | /// \param graph The directed graph the algorithm runs on. |
---|
| 648 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 649 | /// \param cost The cost (length) values of the edges. |
---|
| 650 | /// \param supply The supply values of the nodes (signed). |
---|
| 651 | NetworkSimplex( const Graph &graph, |
---|
| 652 | const CapacityMap &capacity, |
---|
| 653 | const CostMap &cost, |
---|
| 654 | const SupplyMap &supply ) : |
---|
| 655 | _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), |
---|
| 656 | _cost(_graph), _supply(_graph), _flow(_graph), |
---|
| 657 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 658 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 659 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
[2623] | 660 | _flow_result(NULL), _potential_result(NULL), |
---|
[2581] | 661 | _local_flow(false), _local_potential(false), |
---|
[2575] | 662 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 663 | { |
---|
| 664 | // Checking the sum of supply values |
---|
| 665 | Supply sum = 0; |
---|
[2575] | 666 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
---|
| 667 | sum += supply[n]; |
---|
| 668 | if (!(_valid_supply = sum == 0)) return; |
---|
[2440] | 669 | |
---|
[2575] | 670 | // Copying _graph_ref to graph |
---|
| 671 | copyGraph(_graph, _graph_ref) |
---|
| 672 | .edgeMap(_capacity, capacity) |
---|
| 673 | .edgeMap(_cost, cost) |
---|
| 674 | .nodeMap(_supply, supply) |
---|
| 675 | .nodeRef(_node_ref) |
---|
| 676 | .edgeRef(_edge_ref) |
---|
[2556] | 677 | .run(); |
---|
[2440] | 678 | } |
---|
| 679 | |
---|
[2581] | 680 | /// \brief Simple constructor (with lower bounds). |
---|
[2440] | 681 | /// |
---|
[2581] | 682 | /// Simple constructor (with lower bounds). |
---|
[2440] | 683 | /// |
---|
[2575] | 684 | /// \param graph The directed graph the algorithm runs on. |
---|
| 685 | /// \param lower The lower bounds of the edges. |
---|
| 686 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 687 | /// \param cost The cost (length) values of the edges. |
---|
| 688 | /// \param s The source node. |
---|
| 689 | /// \param t The target node. |
---|
| 690 | /// \param flow_value The required amount of flow from node \c s |
---|
| 691 | /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
---|
| 692 | NetworkSimplex( const Graph &graph, |
---|
| 693 | const LowerMap &lower, |
---|
| 694 | const CapacityMap &capacity, |
---|
| 695 | const CostMap &cost, |
---|
| 696 | typename Graph::Node s, |
---|
| 697 | typename Graph::Node t, |
---|
| 698 | typename SupplyMap::Value flow_value ) : |
---|
| 699 | _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), |
---|
| 700 | _cost(_graph), _supply(_graph), _flow(_graph), |
---|
| 701 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 702 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 703 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
[2623] | 704 | _flow_result(NULL), _potential_result(NULL), |
---|
[2581] | 705 | _local_flow(false), _local_potential(false), |
---|
[2575] | 706 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 707 | { |
---|
[2575] | 708 | // Copying _graph_ref to graph |
---|
| 709 | copyGraph(_graph, _graph_ref) |
---|
| 710 | .edgeMap(_cost, cost) |
---|
| 711 | .nodeRef(_node_ref) |
---|
| 712 | .edgeRef(_edge_ref) |
---|
[2556] | 713 | .run(); |
---|
[2440] | 714 | |
---|
[2556] | 715 | // Removing non-zero lower bounds |
---|
[2575] | 716 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { |
---|
| 717 | _capacity[_edge_ref[e]] = capacity[e] - lower[e]; |
---|
[2440] | 718 | } |
---|
[2575] | 719 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { |
---|
| 720 | Supply sum = 0; |
---|
| 721 | if (n == s) sum = flow_value; |
---|
| 722 | if (n == t) sum = -flow_value; |
---|
| 723 | for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 724 | sum += lower[e]; |
---|
| 725 | for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 726 | sum -= lower[e]; |
---|
| 727 | _supply[_node_ref[n]] = sum; |
---|
[2440] | 728 | } |
---|
[2575] | 729 | _valid_supply = true; |
---|
[2440] | 730 | } |
---|
| 731 | |
---|
[2581] | 732 | /// \brief Simple constructor (without lower bounds). |
---|
[2440] | 733 | /// |
---|
[2581] | 734 | /// Simple constructor (without lower bounds). |
---|
[2440] | 735 | /// |
---|
[2575] | 736 | /// \param graph The directed graph the algorithm runs on. |
---|
| 737 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 738 | /// \param cost The cost (length) values of the edges. |
---|
| 739 | /// \param s The source node. |
---|
| 740 | /// \param t The target node. |
---|
| 741 | /// \param flow_value The required amount of flow from node \c s |
---|
| 742 | /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
---|
| 743 | NetworkSimplex( const Graph &graph, |
---|
| 744 | const CapacityMap &capacity, |
---|
| 745 | const CostMap &cost, |
---|
| 746 | typename Graph::Node s, |
---|
| 747 | typename Graph::Node t, |
---|
| 748 | typename SupplyMap::Value flow_value ) : |
---|
| 749 | _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), |
---|
| 750 | _cost(_graph), _supply(_graph, 0), _flow(_graph), |
---|
| 751 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 752 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 753 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
[2623] | 754 | _flow_result(NULL), _potential_result(NULL), |
---|
[2581] | 755 | _local_flow(false), _local_potential(false), |
---|
[2575] | 756 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 757 | { |
---|
[2575] | 758 | // Copying _graph_ref to graph |
---|
| 759 | copyGraph(_graph, _graph_ref) |
---|
| 760 | .edgeMap(_capacity, capacity) |
---|
| 761 | .edgeMap(_cost, cost) |
---|
| 762 | .nodeRef(_node_ref) |
---|
| 763 | .edgeRef(_edge_ref) |
---|
[2556] | 764 | .run(); |
---|
[2575] | 765 | _supply[_node_ref[s]] = flow_value; |
---|
| 766 | _supply[_node_ref[t]] = -flow_value; |
---|
| 767 | _valid_supply = true; |
---|
[2440] | 768 | } |
---|
| 769 | |
---|
[2581] | 770 | /// Destructor. |
---|
| 771 | ~NetworkSimplex() { |
---|
| 772 | if (_local_flow) delete _flow_result; |
---|
| 773 | if (_local_potential) delete _potential_result; |
---|
| 774 | } |
---|
| 775 | |
---|
[2619] | 776 | /// \brief Set the flow map. |
---|
[2581] | 777 | /// |
---|
[2619] | 778 | /// Set the flow map. |
---|
[2581] | 779 | /// |
---|
| 780 | /// \return \c (*this) |
---|
| 781 | NetworkSimplex& flowMap(FlowMap &map) { |
---|
| 782 | if (_local_flow) { |
---|
| 783 | delete _flow_result; |
---|
| 784 | _local_flow = false; |
---|
| 785 | } |
---|
| 786 | _flow_result = ↦ |
---|
| 787 | return *this; |
---|
| 788 | } |
---|
| 789 | |
---|
[2619] | 790 | /// \brief Set the potential map. |
---|
[2581] | 791 | /// |
---|
[2619] | 792 | /// Set the potential map. |
---|
[2581] | 793 | /// |
---|
| 794 | /// \return \c (*this) |
---|
| 795 | NetworkSimplex& potentialMap(PotentialMap &map) { |
---|
| 796 | if (_local_potential) { |
---|
| 797 | delete _potential_result; |
---|
| 798 | _local_potential = false; |
---|
| 799 | } |
---|
| 800 | _potential_result = ↦ |
---|
| 801 | return *this; |
---|
| 802 | } |
---|
| 803 | |
---|
| 804 | /// \name Execution control |
---|
| 805 | |
---|
| 806 | /// @{ |
---|
| 807 | |
---|
[2556] | 808 | /// \brief Runs the algorithm. |
---|
| 809 | /// |
---|
| 810 | /// Runs the algorithm. |
---|
| 811 | /// |
---|
[2575] | 812 | /// \param pivot_rule The pivot rule that is used during the |
---|
| 813 | /// algorithm. |
---|
| 814 | /// |
---|
| 815 | /// The available pivot rules: |
---|
| 816 | /// |
---|
| 817 | /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in |
---|
| 818 | /// a wraparound fashion in every iteration |
---|
| 819 | /// (\ref FirstEligiblePivotRule). |
---|
| 820 | /// |
---|
| 821 | /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in |
---|
| 822 | /// every iteration (\ref BestEligiblePivotRule). |
---|
| 823 | /// |
---|
| 824 | /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in |
---|
| 825 | /// every iteration in a wraparound fashion and the best eligible |
---|
| 826 | /// edge is selected from this block (\ref BlockSearchPivotRule). |
---|
| 827 | /// |
---|
[2619] | 828 | /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is |
---|
| 829 | /// built from eligible edges in a wraparound fashion and in the |
---|
| 830 | /// following minor iterations the best eligible edge is selected |
---|
| 831 | /// from this list (\ref CandidateListPivotRule). |
---|
[2575] | 832 | /// |
---|
[2619] | 833 | /// - ALTERING_LIST_PIVOT It is a modified version of the |
---|
| 834 | /// "Candidate List" pivot rule. It keeps only the several best |
---|
| 835 | /// eligible edges from the former candidate list and extends this |
---|
| 836 | /// list in every iteration (\ref AlteringListPivotRule). |
---|
[2575] | 837 | /// |
---|
[2619] | 838 | /// According to our comprehensive benchmark tests the "Block Search" |
---|
| 839 | /// pivot rule proved to be by far the fastest and the most robust |
---|
| 840 | /// on various test inputs. Thus it is the default option. |
---|
[2575] | 841 | /// |
---|
[2556] | 842 | /// \return \c true if a feasible flow can be found. |
---|
[2619] | 843 | bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) { |
---|
[2575] | 844 | return init() && start(pivot_rule); |
---|
[2556] | 845 | } |
---|
| 846 | |
---|
[2581] | 847 | /// @} |
---|
| 848 | |
---|
| 849 | /// \name Query Functions |
---|
[2619] | 850 | /// The results of the algorithm can be obtained using these |
---|
| 851 | /// functions.\n |
---|
| 852 | /// \ref lemon::NetworkSimplex::run() "run()" must be called before |
---|
| 853 | /// using them. |
---|
[2581] | 854 | |
---|
| 855 | /// @{ |
---|
| 856 | |
---|
[2619] | 857 | /// \brief Return a const reference to the edge map storing the |
---|
[2575] | 858 | /// found flow. |
---|
[2440] | 859 | /// |
---|
[2619] | 860 | /// Return a const reference to the edge map storing the found flow. |
---|
[2440] | 861 | /// |
---|
| 862 | /// \pre \ref run() must be called before using this function. |
---|
| 863 | const FlowMap& flowMap() const { |
---|
[2581] | 864 | return *_flow_result; |
---|
[2440] | 865 | } |
---|
| 866 | |
---|
[2619] | 867 | /// \brief Return a const reference to the node map storing the |
---|
[2575] | 868 | /// found potentials (the dual solution). |
---|
[2440] | 869 | /// |
---|
[2619] | 870 | /// Return a const reference to the node map storing the found |
---|
[2575] | 871 | /// potentials (the dual solution). |
---|
[2440] | 872 | /// |
---|
| 873 | /// \pre \ref run() must be called before using this function. |
---|
| 874 | const PotentialMap& potentialMap() const { |
---|
[2581] | 875 | return *_potential_result; |
---|
| 876 | } |
---|
| 877 | |
---|
[2619] | 878 | /// \brief Return the flow on the given edge. |
---|
[2581] | 879 | /// |
---|
[2619] | 880 | /// Return the flow on the given edge. |
---|
[2581] | 881 | /// |
---|
| 882 | /// \pre \ref run() must be called before using this function. |
---|
| 883 | Capacity flow(const typename Graph::Edge& edge) const { |
---|
| 884 | return (*_flow_result)[edge]; |
---|
| 885 | } |
---|
| 886 | |
---|
[2619] | 887 | /// \brief Return the potential of the given node. |
---|
[2581] | 888 | /// |
---|
[2619] | 889 | /// Return the potential of the given node. |
---|
[2581] | 890 | /// |
---|
| 891 | /// \pre \ref run() must be called before using this function. |
---|
| 892 | Cost potential(const typename Graph::Node& node) const { |
---|
| 893 | return (*_potential_result)[node]; |
---|
[2440] | 894 | } |
---|
| 895 | |
---|
[2619] | 896 | /// \brief Return the total cost of the found flow. |
---|
[2440] | 897 | /// |
---|
[2619] | 898 | /// Return the total cost of the found flow. The complexity of the |
---|
[2440] | 899 | /// function is \f$ O(e) \f$. |
---|
| 900 | /// |
---|
| 901 | /// \pre \ref run() must be called before using this function. |
---|
| 902 | Cost totalCost() const { |
---|
| 903 | Cost c = 0; |
---|
[2575] | 904 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
---|
[2581] | 905 | c += (*_flow_result)[e] * _cost[_edge_ref[e]]; |
---|
[2440] | 906 | return c; |
---|
| 907 | } |
---|
| 908 | |
---|
[2581] | 909 | /// @} |
---|
| 910 | |
---|
[2575] | 911 | private: |
---|
[2440] | 912 | |
---|
[2619] | 913 | /// \brief Extend the underlying graph and initialize all the |
---|
[2440] | 914 | /// node and edge maps. |
---|
| 915 | bool init() { |
---|
[2575] | 916 | if (!_valid_supply) return false; |
---|
[2440] | 917 | |
---|
[2581] | 918 | // Initializing result maps |
---|
| 919 | if (!_flow_result) { |
---|
| 920 | _flow_result = new FlowMap(_graph_ref); |
---|
| 921 | _local_flow = true; |
---|
| 922 | } |
---|
| 923 | if (!_potential_result) { |
---|
| 924 | _potential_result = new PotentialMap(_graph_ref); |
---|
| 925 | _local_potential = true; |
---|
| 926 | } |
---|
| 927 | |
---|
[2619] | 928 | // Initializing the edge vector and the edge maps |
---|
| 929 | _edges.reserve(countEdges(_graph)); |
---|
[2575] | 930 | for (EdgeIt e(_graph); e != INVALID; ++e) { |
---|
[2619] | 931 | _edges.push_back(e); |
---|
[2575] | 932 | _flow[e] = 0; |
---|
| 933 | _state[e] = STATE_LOWER; |
---|
[2440] | 934 | } |
---|
| 935 | |
---|
| 936 | // Adding an artificial root node to the graph |
---|
[2575] | 937 | _root = _graph.addNode(); |
---|
| 938 | _parent[_root] = INVALID; |
---|
| 939 | _pred_edge[_root] = INVALID; |
---|
| 940 | _depth[_root] = 0; |
---|
| 941 | _supply[_root] = 0; |
---|
| 942 | _potential[_root] = 0; |
---|
[2440] | 943 | |
---|
| 944 | // Adding artificial edges to the graph and initializing the node |
---|
| 945 | // maps of the spanning tree data structure |
---|
[2575] | 946 | Node last = _root; |
---|
[2440] | 947 | Edge e; |
---|
| 948 | Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
---|
[2575] | 949 | for (NodeIt u(_graph); u != INVALID; ++u) { |
---|
| 950 | if (u == _root) continue; |
---|
| 951 | _thread[last] = u; |
---|
[2556] | 952 | last = u; |
---|
[2575] | 953 | _parent[u] = _root; |
---|
| 954 | _depth[u] = 1; |
---|
| 955 | if (_supply[u] >= 0) { |
---|
| 956 | e = _graph.addEdge(u, _root); |
---|
| 957 | _flow[e] = _supply[u]; |
---|
| 958 | _forward[u] = true; |
---|
| 959 | _potential[u] = -max_cost; |
---|
[2556] | 960 | } else { |
---|
[2575] | 961 | e = _graph.addEdge(_root, u); |
---|
| 962 | _flow[e] = -_supply[u]; |
---|
| 963 | _forward[u] = false; |
---|
| 964 | _potential[u] = max_cost; |
---|
[2556] | 965 | } |
---|
[2575] | 966 | _cost[e] = max_cost; |
---|
| 967 | _capacity[e] = std::numeric_limits<Capacity>::max(); |
---|
| 968 | _state[e] = STATE_TREE; |
---|
| 969 | _pred_edge[u] = e; |
---|
[2440] | 970 | } |
---|
[2575] | 971 | _thread[last] = _root; |
---|
[2440] | 972 | |
---|
[2575] | 973 | return true; |
---|
[2440] | 974 | } |
---|
| 975 | |
---|
[2619] | 976 | /// Find the join node. |
---|
| 977 | inline Node findJoinNode() { |
---|
[2575] | 978 | Node u = _graph.source(_in_edge); |
---|
| 979 | Node v = _graph.target(_in_edge); |
---|
| 980 | while (u != v) { |
---|
| 981 | if (_depth[u] == _depth[v]) { |
---|
| 982 | u = _parent[u]; |
---|
| 983 | v = _parent[v]; |
---|
[2556] | 984 | } |
---|
[2575] | 985 | else if (_depth[u] > _depth[v]) u = _parent[u]; |
---|
| 986 | else v = _parent[v]; |
---|
[2440] | 987 | } |
---|
| 988 | return u; |
---|
| 989 | } |
---|
| 990 | |
---|
[2619] | 991 | /// \brief Find the leaving edge of the cycle. |
---|
| 992 | /// \return \c true if the leaving edge is not the same as the |
---|
| 993 | /// entering edge. |
---|
| 994 | inline bool findLeavingEdge() { |
---|
[2440] | 995 | // Initializing first and second nodes according to the direction |
---|
| 996 | // of the cycle |
---|
[2575] | 997 | if (_state[_in_edge] == STATE_LOWER) { |
---|
| 998 | first = _graph.source(_in_edge); |
---|
| 999 | second = _graph.target(_in_edge); |
---|
[2440] | 1000 | } else { |
---|
[2575] | 1001 | first = _graph.target(_in_edge); |
---|
| 1002 | second = _graph.source(_in_edge); |
---|
[2440] | 1003 | } |
---|
[2575] | 1004 | delta = _capacity[_in_edge]; |
---|
[2440] | 1005 | bool result = false; |
---|
| 1006 | Capacity d; |
---|
| 1007 | Edge e; |
---|
| 1008 | |
---|
| 1009 | // Searching the cycle along the path form the first node to the |
---|
| 1010 | // root node |
---|
[2575] | 1011 | for (Node u = first; u != join; u = _parent[u]) { |
---|
| 1012 | e = _pred_edge[u]; |
---|
| 1013 | d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e]; |
---|
[2556] | 1014 | if (d < delta) { |
---|
| 1015 | delta = d; |
---|
| 1016 | u_out = u; |
---|
| 1017 | u_in = first; |
---|
| 1018 | v_in = second; |
---|
| 1019 | result = true; |
---|
| 1020 | } |
---|
[2440] | 1021 | } |
---|
| 1022 | // Searching the cycle along the path form the second node to the |
---|
| 1023 | // root node |
---|
[2575] | 1024 | for (Node u = second; u != join; u = _parent[u]) { |
---|
| 1025 | e = _pred_edge[u]; |
---|
| 1026 | d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e]; |
---|
[2556] | 1027 | if (d <= delta) { |
---|
| 1028 | delta = d; |
---|
| 1029 | u_out = u; |
---|
| 1030 | u_in = second; |
---|
| 1031 | v_in = first; |
---|
| 1032 | result = true; |
---|
| 1033 | } |
---|
[2440] | 1034 | } |
---|
| 1035 | return result; |
---|
| 1036 | } |
---|
| 1037 | |
---|
[2619] | 1038 | /// Change \c flow and \c state edge maps. |
---|
| 1039 | inline void changeFlows(bool change) { |
---|
[2440] | 1040 | // Augmenting along the cycle |
---|
| 1041 | if (delta > 0) { |
---|
[2575] | 1042 | Capacity val = _state[_in_edge] * delta; |
---|
| 1043 | _flow[_in_edge] += val; |
---|
| 1044 | for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) { |
---|
| 1045 | _flow[_pred_edge[u]] += _forward[u] ? -val : val; |
---|
[2556] | 1046 | } |
---|
[2575] | 1047 | for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) { |
---|
| 1048 | _flow[_pred_edge[u]] += _forward[u] ? val : -val; |
---|
[2556] | 1049 | } |
---|
[2440] | 1050 | } |
---|
| 1051 | // Updating the state of the entering and leaving edges |
---|
| 1052 | if (change) { |
---|
[2575] | 1053 | _state[_in_edge] = STATE_TREE; |
---|
| 1054 | _state[_pred_edge[u_out]] = |
---|
| 1055 | (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; |
---|
[2440] | 1056 | } else { |
---|
[2575] | 1057 | _state[_in_edge] = -_state[_in_edge]; |
---|
[2440] | 1058 | } |
---|
| 1059 | } |
---|
| 1060 | |
---|
[2619] | 1061 | /// Update \c thread and \c parent node maps. |
---|
| 1062 | inline void updateThreadParent() { |
---|
[2440] | 1063 | Node u; |
---|
[2575] | 1064 | v_out = _parent[u_out]; |
---|
[2440] | 1065 | |
---|
| 1066 | // Handling the case when join and v_out coincide |
---|
| 1067 | bool par_first = false; |
---|
| 1068 | if (join == v_out) { |
---|
[2575] | 1069 | for (u = join; u != u_in && u != v_in; u = _thread[u]) ; |
---|
[2556] | 1070 | if (u == v_in) { |
---|
| 1071 | par_first = true; |
---|
[2575] | 1072 | while (_thread[u] != u_out) u = _thread[u]; |
---|
[2556] | 1073 | first = u; |
---|
| 1074 | } |
---|
[2440] | 1075 | } |
---|
| 1076 | |
---|
| 1077 | // Finding the last successor of u_in (u) and the node after it |
---|
| 1078 | // (right) according to the thread index |
---|
[2575] | 1079 | for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ; |
---|
| 1080 | right = _thread[u]; |
---|
| 1081 | if (_thread[v_in] == u_out) { |
---|
| 1082 | for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ; |
---|
| 1083 | if (last == u_out) last = _thread[last]; |
---|
[2440] | 1084 | } |
---|
[2575] | 1085 | else last = _thread[v_in]; |
---|
[2440] | 1086 | |
---|
| 1087 | // Updating stem nodes |
---|
[2575] | 1088 | _thread[v_in] = stem = u_in; |
---|
[2440] | 1089 | par_stem = v_in; |
---|
| 1090 | while (stem != u_out) { |
---|
[2575] | 1091 | _thread[u] = new_stem = _parent[stem]; |
---|
[2440] | 1092 | |
---|
[2556] | 1093 | // Finding the node just before the stem node (u) according to |
---|
| 1094 | // the original thread index |
---|
[2575] | 1095 | for (u = new_stem; _thread[u] != stem; u = _thread[u]) ; |
---|
| 1096 | _thread[u] = right; |
---|
[2440] | 1097 | |
---|
[2556] | 1098 | // Changing the parent node of stem and shifting stem and |
---|
| 1099 | // par_stem nodes |
---|
[2575] | 1100 | _parent[stem] = par_stem; |
---|
[2556] | 1101 | par_stem = stem; |
---|
| 1102 | stem = new_stem; |
---|
[2440] | 1103 | |
---|
[2556] | 1104 | // Finding the last successor of stem (u) and the node after it |
---|
| 1105 | // (right) according to the thread index |
---|
[2575] | 1106 | for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ; |
---|
| 1107 | right = _thread[u]; |
---|
[2440] | 1108 | } |
---|
[2575] | 1109 | _parent[u_out] = par_stem; |
---|
| 1110 | _thread[u] = last; |
---|
[2440] | 1111 | |
---|
| 1112 | if (join == v_out && par_first) { |
---|
[2575] | 1113 | if (first != v_in) _thread[first] = right; |
---|
[2440] | 1114 | } else { |
---|
[2575] | 1115 | for (u = v_out; _thread[u] != u_out; u = _thread[u]) ; |
---|
| 1116 | _thread[u] = right; |
---|
[2440] | 1117 | } |
---|
| 1118 | } |
---|
| 1119 | |
---|
[2619] | 1120 | /// Update \c pred_edge and \c forward node maps. |
---|
| 1121 | inline void updatePredEdge() { |
---|
[2440] | 1122 | Node u = u_out, v; |
---|
| 1123 | while (u != u_in) { |
---|
[2575] | 1124 | v = _parent[u]; |
---|
| 1125 | _pred_edge[u] = _pred_edge[v]; |
---|
| 1126 | _forward[u] = !_forward[v]; |
---|
[2556] | 1127 | u = v; |
---|
[2440] | 1128 | } |
---|
[2575] | 1129 | _pred_edge[u_in] = _in_edge; |
---|
| 1130 | _forward[u_in] = (u_in == _graph.source(_in_edge)); |
---|
[2440] | 1131 | } |
---|
| 1132 | |
---|
[2619] | 1133 | /// Update \c depth and \c potential node maps. |
---|
| 1134 | inline void updateDepthPotential() { |
---|
[2575] | 1135 | _depth[u_in] = _depth[v_in] + 1; |
---|
| 1136 | _potential[u_in] = _forward[u_in] ? |
---|
| 1137 | _potential[v_in] - _cost[_pred_edge[u_in]] : |
---|
| 1138 | _potential[v_in] + _cost[_pred_edge[u_in]]; |
---|
[2440] | 1139 | |
---|
[2575] | 1140 | Node u = _thread[u_in], v; |
---|
[2440] | 1141 | while (true) { |
---|
[2575] | 1142 | v = _parent[u]; |
---|
[2556] | 1143 | if (v == INVALID) break; |
---|
[2575] | 1144 | _depth[u] = _depth[v] + 1; |
---|
| 1145 | _potential[u] = _forward[u] ? |
---|
| 1146 | _potential[v] - _cost[_pred_edge[u]] : |
---|
| 1147 | _potential[v] + _cost[_pred_edge[u]]; |
---|
| 1148 | if (_depth[u] <= _depth[v_in]) break; |
---|
| 1149 | u = _thread[u]; |
---|
[2440] | 1150 | } |
---|
| 1151 | } |
---|
| 1152 | |
---|
[2619] | 1153 | /// Execute the algorithm. |
---|
[2575] | 1154 | bool start(PivotRuleEnum pivot_rule) { |
---|
[2619] | 1155 | // Selecting the pivot rule implementation |
---|
[2575] | 1156 | switch (pivot_rule) { |
---|
| 1157 | case FIRST_ELIGIBLE_PIVOT: |
---|
| 1158 | return start<FirstEligiblePivotRule>(); |
---|
| 1159 | case BEST_ELIGIBLE_PIVOT: |
---|
| 1160 | return start<BestEligiblePivotRule>(); |
---|
| 1161 | case BLOCK_SEARCH_PIVOT: |
---|
| 1162 | return start<BlockSearchPivotRule>(); |
---|
| 1163 | case CANDIDATE_LIST_PIVOT: |
---|
| 1164 | return start<CandidateListPivotRule>(); |
---|
[2619] | 1165 | case ALTERING_LIST_PIVOT: |
---|
| 1166 | return start<AlteringListPivotRule>(); |
---|
[2575] | 1167 | } |
---|
| 1168 | return false; |
---|
| 1169 | } |
---|
| 1170 | |
---|
| 1171 | template<class PivotRuleImplementation> |
---|
[2440] | 1172 | bool start() { |
---|
[2619] | 1173 | PivotRuleImplementation pivot(*this, _edges); |
---|
[2575] | 1174 | |
---|
| 1175 | // Executing the network simplex algorithm |
---|
| 1176 | while (pivot.findEnteringEdge()) { |
---|
[2556] | 1177 | join = findJoinNode(); |
---|
| 1178 | bool change = findLeavingEdge(); |
---|
| 1179 | changeFlows(change); |
---|
| 1180 | if (change) { |
---|
| 1181 | updateThreadParent(); |
---|
| 1182 | updatePredEdge(); |
---|
| 1183 | updateDepthPotential(); |
---|
| 1184 | } |
---|
[2440] | 1185 | } |
---|
| 1186 | |
---|
[2575] | 1187 | // Checking if the flow amount equals zero on all the artificial |
---|
| 1188 | // edges |
---|
| 1189 | for (InEdgeIt e(_graph, _root); e != INVALID; ++e) |
---|
| 1190 | if (_flow[e] > 0) return false; |
---|
| 1191 | for (OutEdgeIt e(_graph, _root); e != INVALID; ++e) |
---|
| 1192 | if (_flow[e] > 0) return false; |
---|
[2440] | 1193 | |
---|
[2575] | 1194 | // Copying flow values to _flow_result |
---|
| 1195 | if (_lower) { |
---|
| 1196 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
---|
[2581] | 1197 | (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]]; |
---|
[2440] | 1198 | } else { |
---|
[2575] | 1199 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
---|
[2581] | 1200 | (*_flow_result)[e] = _flow[_edge_ref[e]]; |
---|
[2440] | 1201 | } |
---|
[2575] | 1202 | // Copying potential values to _potential_result |
---|
| 1203 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
---|
[2581] | 1204 | (*_potential_result)[n] = _potential[_node_ref[n]]; |
---|
[2440] | 1205 | |
---|
| 1206 | return true; |
---|
| 1207 | } |
---|
| 1208 | |
---|
| 1209 | }; //class NetworkSimplex |
---|
| 1210 | |
---|
| 1211 | ///@} |
---|
| 1212 | |
---|
| 1213 | } //namespace lemon |
---|
| 1214 | |
---|
| 1215 | #endif //LEMON_NETWORK_SIMPLEX_H |
---|