deba@2440: /* -*- C++ -*- deba@2440: * deba@2440: * This file is a part of LEMON, a generic C++ optimization library deba@2440: * alpar@2553: * Copyright (C) 2003-2008 deba@2440: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport deba@2440: * (Egervary Research Group on Combinatorial Optimization, EGRES). deba@2440: * deba@2440: * Permission to use, modify and distribute this software is granted deba@2440: * provided that this copyright notice appears in all copies. For deba@2440: * precise terms see the accompanying LICENSE file. deba@2440: * deba@2440: * This software is provided "AS IS" with no warranty of any kind, deba@2440: * express or implied, and with no claim as to its suitability for any deba@2440: * purpose. deba@2440: * deba@2440: */ deba@2440: deba@2440: #ifndef LEMON_NETWORK_SIMPLEX_H deba@2440: #define LEMON_NETWORK_SIMPLEX_H deba@2440: deba@2440: /// \ingroup min_cost_flow deba@2440: /// deba@2440: /// \file kpeter@2575: /// \brief Network simplex algorithm for finding a minimum cost flow. deba@2440: kpeter@2575: #include deba@2440: #include kpeter@2575: kpeter@2509: #include kpeter@2509: #include deba@2440: #include kpeter@2575: #include deba@2440: deba@2440: namespace lemon { deba@2440: deba@2440: /// \addtogroup min_cost_flow deba@2440: /// @{ deba@2440: deba@2440: /// \brief Implementation of the network simplex algorithm for deba@2440: /// finding a minimum cost flow. deba@2440: /// kpeter@2556: /// \ref NetworkSimplex implements the network simplex algorithm for kpeter@2556: /// finding a minimum cost flow. deba@2440: /// kpeter@2575: /// \tparam Graph The directed graph type the algorithm runs on. kpeter@2575: /// \tparam LowerMap The type of the lower bound map. kpeter@2575: /// \tparam CapacityMap The type of the capacity (upper bound) map. kpeter@2575: /// \tparam CostMap The type of the cost (length) map. kpeter@2575: /// \tparam SupplyMap The type of the supply map. deba@2440: /// deba@2440: /// \warning kpeter@2575: /// - Edge capacities and costs should be \e non-negative \e integers. kpeter@2575: /// - Supply values should be \e signed \e integers. kpeter@2581: /// - The value types of the maps should be convertible to each other. kpeter@2581: /// - \c CostMap::Value must be signed type. kpeter@2575: /// kpeter@2575: /// \note \ref NetworkSimplex provides six different pivot rule kpeter@2575: /// implementations that significantly affect the efficiency of the kpeter@2575: /// algorithm. kpeter@2575: /// By default a combined pivot rule is used, which is the fastest kpeter@2575: /// implementation according to our benchmark tests. kpeter@2575: /// Another pivot rule can be selected using \ref run() function kpeter@2575: /// with the proper parameter. deba@2440: /// deba@2440: /// \author Peter Kovacs deba@2440: kpeter@2533: template < typename Graph, kpeter@2533: typename LowerMap = typename Graph::template EdgeMap, kpeter@2575: typename CapacityMap = typename Graph::template EdgeMap, kpeter@2533: typename CostMap = typename Graph::template EdgeMap, kpeter@2575: typename SupplyMap = typename Graph::template NodeMap > deba@2440: class NetworkSimplex deba@2440: { deba@2440: typedef typename CapacityMap::Value Capacity; deba@2440: typedef typename CostMap::Value Cost; deba@2440: typedef typename SupplyMap::Value Supply; deba@2440: deba@2440: typedef SmartGraph SGraph; kpeter@2556: GRAPH_TYPEDEFS(typename SGraph); deba@2440: deba@2440: typedef typename SGraph::template EdgeMap SCapacityMap; deba@2440: typedef typename SGraph::template EdgeMap SCostMap; deba@2440: typedef typename SGraph::template NodeMap SSupplyMap; deba@2440: typedef typename SGraph::template NodeMap SPotentialMap; deba@2440: deba@2440: typedef typename SGraph::template NodeMap IntNodeMap; deba@2440: typedef typename SGraph::template NodeMap BoolNodeMap; deba@2440: typedef typename SGraph::template NodeMap NodeNodeMap; deba@2440: typedef typename SGraph::template NodeMap EdgeNodeMap; deba@2440: typedef typename SGraph::template EdgeMap IntEdgeMap; deba@2440: deba@2440: typedef typename Graph::template NodeMap NodeRefMap; deba@2440: typedef typename Graph::template EdgeMap EdgeRefMap; deba@2440: deba@2440: public: deba@2440: kpeter@2556: /// The type of the flow map. deba@2440: typedef typename Graph::template EdgeMap FlowMap; kpeter@2556: /// The type of the potential map. deba@2440: typedef typename Graph::template NodeMap PotentialMap; deba@2440: kpeter@2575: public: deba@2440: kpeter@2575: /// Enum type to select the pivot rule used by \ref run(). kpeter@2575: enum PivotRuleEnum { kpeter@2575: FIRST_ELIGIBLE_PIVOT, kpeter@2575: BEST_ELIGIBLE_PIVOT, kpeter@2575: BLOCK_SEARCH_PIVOT, kpeter@2575: LIMITED_SEARCH_PIVOT, kpeter@2575: CANDIDATE_LIST_PIVOT, kpeter@2575: COMBINED_PIVOT kpeter@2575: }; kpeter@2575: kpeter@2575: private: kpeter@2575: kpeter@2575: /// \brief Map adaptor class for handling reduced edge costs. kpeter@2575: /// kpeter@2556: /// Map adaptor class for handling reduced edge costs. deba@2440: class ReducedCostMap : public MapBase deba@2440: { deba@2440: private: deba@2440: kpeter@2575: const SGraph &_gr; kpeter@2575: const SCostMap &_cost_map; kpeter@2575: const SPotentialMap &_pot_map; deba@2440: deba@2440: public: deba@2440: kpeter@2575: ///\e kpeter@2575: ReducedCostMap( const SGraph &gr, kpeter@2575: const SCostMap &cost_map, kpeter@2575: const SPotentialMap &pot_map ) : kpeter@2579: _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {} deba@2440: kpeter@2575: ///\e kpeter@2509: Cost operator[](const Edge &e) const { kpeter@2575: return _cost_map[e] + _pot_map[_gr.source(e)] kpeter@2575: - _pot_map[_gr.target(e)]; deba@2440: } deba@2440: deba@2440: }; //class ReducedCostMap deba@2440: kpeter@2575: private: deba@2440: kpeter@2575: /// \brief Implementation of the "First Eligible" pivot rule for the kpeter@2575: /// \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: /// kpeter@2575: /// This class implements the "First Eligible" pivot rule kpeter@2575: /// for the \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: class FirstEligiblePivotRule kpeter@2575: { kpeter@2575: private: deba@2440: kpeter@2575: NetworkSimplex &_ns; kpeter@2575: EdgeIt _next_edge; deba@2440: kpeter@2575: public: deba@2440: kpeter@2575: /// Constructor. kpeter@2575: FirstEligiblePivotRule(NetworkSimplex &ns) : kpeter@2575: _ns(ns), _next_edge(ns._graph) {} kpeter@2575: kpeter@2575: /// Finds the next entering edge. kpeter@2575: bool findEnteringEdge() { kpeter@2575: for (EdgeIt e = _next_edge; e != INVALID; ++e) { kpeter@2575: if (_ns._state[e] * _ns._red_cost[e] < 0) { kpeter@2575: _ns._in_edge = e; kpeter@2575: _next_edge = ++e; kpeter@2575: return true; kpeter@2575: } kpeter@2575: } kpeter@2575: for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { kpeter@2575: if (_ns._state[e] * _ns._red_cost[e] < 0) { kpeter@2575: _ns._in_edge = e; kpeter@2575: _next_edge = ++e; kpeter@2575: return true; kpeter@2575: } kpeter@2575: } kpeter@2575: return false; kpeter@2575: } kpeter@2575: }; //class FirstEligiblePivotRule kpeter@2575: kpeter@2575: /// \brief Implementation of the "Best Eligible" pivot rule for the kpeter@2575: /// \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: /// kpeter@2575: /// This class implements the "Best Eligible" pivot rule kpeter@2575: /// for the \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: class BestEligiblePivotRule kpeter@2575: { kpeter@2575: private: kpeter@2575: kpeter@2575: NetworkSimplex &_ns; kpeter@2575: kpeter@2575: public: kpeter@2575: kpeter@2575: /// Constructor. kpeter@2575: BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {} kpeter@2575: kpeter@2575: /// Finds the next entering edge. kpeter@2575: bool findEnteringEdge() { kpeter@2575: Cost min = 0; kpeter@2575: for (EdgeIt e(_ns._graph); e != INVALID; ++e) { kpeter@2575: if (_ns._state[e] * _ns._red_cost[e] < min) { kpeter@2575: min = _ns._state[e] * _ns._red_cost[e]; kpeter@2575: _ns._in_edge = e; kpeter@2575: } kpeter@2575: } kpeter@2575: return min < 0; kpeter@2575: } kpeter@2575: }; //class BestEligiblePivotRule kpeter@2575: kpeter@2575: /// \brief Implementation of the "Block Search" pivot rule for the kpeter@2575: /// \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: /// kpeter@2575: /// This class implements the "Block Search" pivot rule kpeter@2575: /// for the \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: class BlockSearchPivotRule kpeter@2575: { kpeter@2575: private: kpeter@2575: kpeter@2575: NetworkSimplex &_ns; kpeter@2575: EdgeIt _next_edge, _min_edge; kpeter@2575: int _block_size; kpeter@2575: kpeter@2575: static const int MIN_BLOCK_SIZE = 10; kpeter@2575: kpeter@2575: public: kpeter@2575: kpeter@2575: /// Constructor. kpeter@2575: BlockSearchPivotRule(NetworkSimplex &ns) : kpeter@2575: _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) kpeter@2575: { kpeter@2575: _block_size = 2 * int(sqrt(countEdges(_ns._graph))); kpeter@2575: if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE; kpeter@2575: } kpeter@2575: kpeter@2575: /// Finds the next entering edge. kpeter@2575: bool findEnteringEdge() { kpeter@2575: Cost curr, min = 0; kpeter@2575: int cnt = 0; kpeter@2575: for (EdgeIt e = _next_edge; e != INVALID; ++e) { kpeter@2575: if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { kpeter@2575: min = curr; kpeter@2575: _min_edge = e; kpeter@2575: } kpeter@2575: if (++cnt == _block_size) { kpeter@2575: if (min < 0) break; kpeter@2575: cnt = 0; kpeter@2575: } kpeter@2575: } kpeter@2575: if (min == 0) { kpeter@2575: for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { kpeter@2575: if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { kpeter@2575: min = curr; kpeter@2575: _min_edge = e; kpeter@2575: } kpeter@2575: if (++cnt == _block_size) { kpeter@2575: if (min < 0) break; kpeter@2575: cnt = 0; kpeter@2575: } kpeter@2575: } kpeter@2575: } kpeter@2575: _ns._in_edge = _min_edge; kpeter@2575: _next_edge = ++_min_edge; kpeter@2575: return min < 0; kpeter@2575: } kpeter@2575: }; //class BlockSearchPivotRule kpeter@2575: kpeter@2575: /// \brief Implementation of the "Limited Search" pivot rule for the kpeter@2575: /// \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: /// kpeter@2575: /// This class implements the "Limited Search" pivot rule kpeter@2575: /// for the \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: class LimitedSearchPivotRule kpeter@2575: { kpeter@2575: private: kpeter@2575: kpeter@2575: NetworkSimplex &_ns; kpeter@2575: EdgeIt _next_edge, _min_edge; kpeter@2575: int _sample_size; kpeter@2575: kpeter@2575: static const int MIN_SAMPLE_SIZE = 10; kpeter@2575: static const double SAMPLE_SIZE_FACTOR = 0.0015; kpeter@2575: kpeter@2575: public: kpeter@2575: kpeter@2575: /// Constructor. kpeter@2575: LimitedSearchPivotRule(NetworkSimplex &ns) : kpeter@2575: _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) kpeter@2575: { kpeter@2575: _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph)); kpeter@2575: if (_sample_size < MIN_SAMPLE_SIZE) kpeter@2575: _sample_size = MIN_SAMPLE_SIZE; kpeter@2575: } kpeter@2575: kpeter@2575: /// Finds the next entering edge. kpeter@2575: bool findEnteringEdge() { kpeter@2575: Cost curr, min = 0; kpeter@2575: int cnt = 0; kpeter@2575: for (EdgeIt e = _next_edge; e != INVALID; ++e) { kpeter@2575: if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { kpeter@2575: min = curr; kpeter@2575: _min_edge = e; kpeter@2575: } kpeter@2575: if (curr < 0 && ++cnt == _sample_size) break; kpeter@2575: } kpeter@2575: if (min == 0) { kpeter@2575: for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { kpeter@2575: if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { kpeter@2575: min = curr; kpeter@2575: _min_edge = e; kpeter@2575: } kpeter@2575: if (curr < 0 && ++cnt == _sample_size) break; kpeter@2575: } kpeter@2575: } kpeter@2575: _ns._in_edge = _min_edge; kpeter@2575: _next_edge = ++_min_edge; kpeter@2575: return min < 0; kpeter@2575: } kpeter@2575: }; //class LimitedSearchPivotRule kpeter@2575: kpeter@2575: /// \brief Implementation of the "Candidate List" pivot rule for the kpeter@2575: /// \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: /// kpeter@2575: /// This class implements the "Candidate List" pivot rule kpeter@2575: /// for the \ref NetworkSimplex "network simplex" algorithm. kpeter@2575: class CandidateListPivotRule kpeter@2575: { kpeter@2575: private: kpeter@2575: kpeter@2575: NetworkSimplex &_ns; kpeter@2575: kpeter@2575: // The list of candidate edges. kpeter@2575: std::vector _candidates; kpeter@2575: // The maximum length of the edge list. kpeter@2575: int _list_length; kpeter@2575: // The maximum number of minor iterations between two major kpeter@2575: // itarations. kpeter@2575: int _minor_limit; kpeter@2575: kpeter@2575: int _minor_count; kpeter@2575: EdgeIt _next_edge; kpeter@2575: kpeter@2575: static const double LIST_LENGTH_FACTOR = 0.002; kpeter@2575: static const double MINOR_LIMIT_FACTOR = 0.1; kpeter@2575: static const int MIN_LIST_LENGTH = 10; kpeter@2575: static const int MIN_MINOR_LIMIT = 2; kpeter@2575: kpeter@2575: public: kpeter@2575: kpeter@2575: /// Constructor. kpeter@2575: CandidateListPivotRule(NetworkSimplex &ns) : kpeter@2575: _ns(ns), _next_edge(ns._graph) kpeter@2575: { kpeter@2575: int edge_num = countEdges(_ns._graph); kpeter@2575: _minor_count = 0; kpeter@2575: _list_length = int(edge_num * LIST_LENGTH_FACTOR); kpeter@2575: if (_list_length < MIN_LIST_LENGTH) kpeter@2575: _list_length = MIN_LIST_LENGTH; kpeter@2575: _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR); kpeter@2575: if (_minor_limit < MIN_MINOR_LIMIT) kpeter@2575: _minor_limit = MIN_MINOR_LIMIT; kpeter@2575: } kpeter@2575: kpeter@2575: /// Finds the next entering edge. kpeter@2575: bool findEnteringEdge() { kpeter@2575: Cost min, curr; kpeter@2575: if (_minor_count < _minor_limit && _candidates.size() > 0) { kpeter@2575: // Minor iteration kpeter@2575: ++_minor_count; kpeter@2575: Edge e; kpeter@2575: min = 0; kpeter@2575: for (int i = 0; i < int(_candidates.size()); ++i) { kpeter@2575: e = _candidates[i]; kpeter@2575: if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { kpeter@2575: min = curr; kpeter@2575: _ns._in_edge = e; kpeter@2575: } kpeter@2575: } kpeter@2575: if (min < 0) return true; kpeter@2575: } kpeter@2575: kpeter@2575: // Major iteration kpeter@2575: _candidates.clear(); kpeter@2575: EdgeIt e = _next_edge; kpeter@2575: min = 0; kpeter@2575: for ( ; e != INVALID; ++e) { kpeter@2575: if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { kpeter@2575: _candidates.push_back(e); kpeter@2575: if (curr < min) { kpeter@2575: min = curr; kpeter@2575: _ns._in_edge = e; kpeter@2575: } kpeter@2575: if (int(_candidates.size()) == _list_length) break; kpeter@2575: } kpeter@2575: } kpeter@2575: if (int(_candidates.size()) < _list_length) { kpeter@2575: for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) { kpeter@2575: if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { kpeter@2575: _candidates.push_back(e); kpeter@2575: if (curr < min) { kpeter@2575: min = curr; kpeter@2575: _ns._in_edge = e; kpeter@2575: } kpeter@2575: if (int(_candidates.size()) == _list_length) break; kpeter@2575: } kpeter@2575: } kpeter@2575: } kpeter@2575: if (_candidates.size() == 0) return false; kpeter@2575: _minor_count = 1; kpeter@2575: _next_edge = ++e; kpeter@2575: return true; kpeter@2575: } kpeter@2575: }; //class CandidateListPivotRule kpeter@2575: kpeter@2575: private: kpeter@2575: kpeter@2579: // State constants for edges kpeter@2579: enum EdgeStateEnum { kpeter@2579: STATE_UPPER = -1, kpeter@2579: STATE_TREE = 0, kpeter@2579: STATE_LOWER = 1 kpeter@2579: }; kpeter@2575: kpeter@2575: // Constant for the combined pivot rule. kpeter@2575: static const int COMBINED_PIVOT_MAX_DEG = 5; kpeter@2575: kpeter@2575: private: kpeter@2575: kpeter@2575: // The directed graph the algorithm runs on kpeter@2575: SGraph _graph; kpeter@2575: // The original graph kpeter@2575: const Graph &_graph_ref; kpeter@2575: // The original lower bound map kpeter@2575: const LowerMap *_lower; kpeter@2575: // The capacity map kpeter@2575: SCapacityMap _capacity; kpeter@2575: // The cost map kpeter@2575: SCostMap _cost; kpeter@2575: // The supply map kpeter@2575: SSupplyMap _supply; kpeter@2575: bool _valid_supply; kpeter@2575: kpeter@2575: // Edge map of the current flow kpeter@2575: SCapacityMap _flow; kpeter@2575: // Node map of the current potentials kpeter@2575: SPotentialMap _potential; kpeter@2575: kpeter@2575: // The depth node map of the spanning tree structure kpeter@2575: IntNodeMap _depth; kpeter@2575: // The parent node map of the spanning tree structure kpeter@2575: NodeNodeMap _parent; kpeter@2575: // The pred_edge node map of the spanning tree structure kpeter@2575: EdgeNodeMap _pred_edge; kpeter@2575: // The thread node map of the spanning tree structure kpeter@2575: NodeNodeMap _thread; kpeter@2575: // The forward node map of the spanning tree structure kpeter@2575: BoolNodeMap _forward; kpeter@2575: // The state edge map kpeter@2575: IntEdgeMap _state; kpeter@2575: // The root node of the starting spanning tree kpeter@2575: Node _root; kpeter@2575: kpeter@2575: // The reduced cost map kpeter@2575: ReducedCostMap _red_cost; kpeter@2575: kpeter@2575: // Members for handling the original graph kpeter@2581: FlowMap *_flow_result; kpeter@2581: PotentialMap *_potential_result; kpeter@2581: bool _local_flow; kpeter@2581: bool _local_potential; kpeter@2575: NodeRefMap _node_ref; kpeter@2575: EdgeRefMap _edge_ref; deba@2440: kpeter@2556: // The entering edge of the current pivot iteration. kpeter@2575: Edge _in_edge; kpeter@2575: kpeter@2556: // Temporary nodes used in the current pivot iteration. kpeter@2556: Node join, u_in, v_in, u_out, v_out; kpeter@2556: Node right, first, second, last; kpeter@2556: Node stem, par_stem, new_stem; kpeter@2556: // The maximum augment amount along the found cycle in the current kpeter@2556: // pivot iteration. kpeter@2556: Capacity delta; deba@2440: deba@2440: public : deba@2440: kpeter@2581: /// \brief General constructor (with lower bounds). deba@2440: /// kpeter@2581: /// General constructor (with lower bounds). deba@2440: /// kpeter@2575: /// \param graph The directed graph the algorithm runs on. kpeter@2575: /// \param lower The lower bounds of the edges. kpeter@2575: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2575: /// \param cost The cost (length) values of the edges. kpeter@2575: /// \param supply The supply values of the nodes (signed). kpeter@2575: NetworkSimplex( const Graph &graph, kpeter@2575: const LowerMap &lower, kpeter@2575: const CapacityMap &capacity, kpeter@2575: const CostMap &cost, kpeter@2575: const SupplyMap &supply ) : kpeter@2575: _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), kpeter@2575: _cost(_graph), _supply(_graph), _flow(_graph), kpeter@2575: _potential(_graph), _depth(_graph), _parent(_graph), kpeter@2575: _pred_edge(_graph), _thread(_graph), _forward(_graph), kpeter@2575: _state(_graph), _red_cost(_graph, _cost, _potential), kpeter@2581: _flow_result(0), _potential_result(0), kpeter@2581: _local_flow(false), _local_potential(false), kpeter@2575: _node_ref(graph), _edge_ref(graph) deba@2440: { deba@2440: // Checking the sum of supply values deba@2440: Supply sum = 0; kpeter@2575: for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) kpeter@2575: sum += supply[n]; kpeter@2575: if (!(_valid_supply = sum == 0)) return; deba@2440: kpeter@2575: // Copying _graph_ref to _graph kpeter@2575: _graph.reserveNode(countNodes(_graph_ref) + 1); kpeter@2575: _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); kpeter@2575: copyGraph(_graph, _graph_ref) kpeter@2575: .edgeMap(_cost, cost) kpeter@2575: .nodeRef(_node_ref) kpeter@2575: .edgeRef(_edge_ref) kpeter@2556: .run(); deba@2440: kpeter@2556: // Removing non-zero lower bounds kpeter@2575: for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { kpeter@2575: _capacity[_edge_ref[e]] = capacity[e] - lower[e]; deba@2440: } kpeter@2575: for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { kpeter@2575: Supply s = supply[n]; kpeter@2575: for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) kpeter@2575: s += lower[e]; kpeter@2575: for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) kpeter@2575: s -= lower[e]; kpeter@2575: _supply[_node_ref[n]] = s; deba@2440: } deba@2440: } deba@2440: kpeter@2581: /// \brief General constructor (without lower bounds). deba@2440: /// kpeter@2581: /// General constructor (without lower bounds). deba@2440: /// kpeter@2575: /// \param graph The directed graph the algorithm runs on. kpeter@2575: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2575: /// \param cost The cost (length) values of the edges. kpeter@2575: /// \param supply The supply values of the nodes (signed). kpeter@2575: NetworkSimplex( const Graph &graph, kpeter@2575: const CapacityMap &capacity, kpeter@2575: const CostMap &cost, kpeter@2575: const SupplyMap &supply ) : kpeter@2575: _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), kpeter@2575: _cost(_graph), _supply(_graph), _flow(_graph), kpeter@2575: _potential(_graph), _depth(_graph), _parent(_graph), kpeter@2575: _pred_edge(_graph), _thread(_graph), _forward(_graph), kpeter@2575: _state(_graph), _red_cost(_graph, _cost, _potential), kpeter@2581: _flow_result(0), _potential_result(0), kpeter@2581: _local_flow(false), _local_potential(false), kpeter@2575: _node_ref(graph), _edge_ref(graph) deba@2440: { deba@2440: // Checking the sum of supply values deba@2440: Supply sum = 0; kpeter@2575: for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) kpeter@2575: sum += supply[n]; kpeter@2575: if (!(_valid_supply = sum == 0)) return; deba@2440: kpeter@2575: // Copying _graph_ref to graph kpeter@2575: copyGraph(_graph, _graph_ref) kpeter@2575: .edgeMap(_capacity, capacity) kpeter@2575: .edgeMap(_cost, cost) kpeter@2575: .nodeMap(_supply, supply) kpeter@2575: .nodeRef(_node_ref) kpeter@2575: .edgeRef(_edge_ref) kpeter@2556: .run(); deba@2440: } deba@2440: kpeter@2581: /// \brief Simple constructor (with lower bounds). deba@2440: /// kpeter@2581: /// Simple constructor (with lower bounds). deba@2440: /// kpeter@2575: /// \param graph The directed graph the algorithm runs on. kpeter@2575: /// \param lower The lower bounds of the edges. kpeter@2575: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2575: /// \param cost The cost (length) values of the edges. kpeter@2575: /// \param s The source node. kpeter@2575: /// \param t The target node. kpeter@2575: /// \param flow_value The required amount of flow from node \c s kpeter@2575: /// to node \c t (i.e. the supply of \c s and the demand of \c t). kpeter@2575: NetworkSimplex( const Graph &graph, kpeter@2575: const LowerMap &lower, kpeter@2575: const CapacityMap &capacity, kpeter@2575: const CostMap &cost, kpeter@2575: typename Graph::Node s, kpeter@2575: typename Graph::Node t, kpeter@2575: typename SupplyMap::Value flow_value ) : kpeter@2575: _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), kpeter@2575: _cost(_graph), _supply(_graph), _flow(_graph), kpeter@2575: _potential(_graph), _depth(_graph), _parent(_graph), kpeter@2575: _pred_edge(_graph), _thread(_graph), _forward(_graph), kpeter@2575: _state(_graph), _red_cost(_graph, _cost, _potential), kpeter@2581: _flow_result(0), _potential_result(0), kpeter@2581: _local_flow(false), _local_potential(false), kpeter@2575: _node_ref(graph), _edge_ref(graph) deba@2440: { kpeter@2575: // Copying _graph_ref to graph kpeter@2575: copyGraph(_graph, _graph_ref) kpeter@2575: .edgeMap(_cost, cost) kpeter@2575: .nodeRef(_node_ref) kpeter@2575: .edgeRef(_edge_ref) kpeter@2556: .run(); deba@2440: kpeter@2556: // Removing non-zero lower bounds kpeter@2575: for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { kpeter@2575: _capacity[_edge_ref[e]] = capacity[e] - lower[e]; deba@2440: } kpeter@2575: for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { kpeter@2575: Supply sum = 0; kpeter@2575: if (n == s) sum = flow_value; kpeter@2575: if (n == t) sum = -flow_value; kpeter@2575: for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) kpeter@2575: sum += lower[e]; kpeter@2575: for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) kpeter@2575: sum -= lower[e]; kpeter@2575: _supply[_node_ref[n]] = sum; deba@2440: } kpeter@2575: _valid_supply = true; deba@2440: } deba@2440: kpeter@2581: /// \brief Simple constructor (without lower bounds). deba@2440: /// kpeter@2581: /// Simple constructor (without lower bounds). deba@2440: /// kpeter@2575: /// \param graph The directed graph the algorithm runs on. kpeter@2575: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2575: /// \param cost The cost (length) values of the edges. kpeter@2575: /// \param s The source node. kpeter@2575: /// \param t The target node. kpeter@2575: /// \param flow_value The required amount of flow from node \c s kpeter@2575: /// to node \c t (i.e. the supply of \c s and the demand of \c t). kpeter@2575: NetworkSimplex( const Graph &graph, kpeter@2575: const CapacityMap &capacity, kpeter@2575: const CostMap &cost, kpeter@2575: typename Graph::Node s, kpeter@2575: typename Graph::Node t, kpeter@2575: typename SupplyMap::Value flow_value ) : kpeter@2575: _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), kpeter@2575: _cost(_graph), _supply(_graph, 0), _flow(_graph), kpeter@2575: _potential(_graph), _depth(_graph), _parent(_graph), kpeter@2575: _pred_edge(_graph), _thread(_graph), _forward(_graph), kpeter@2575: _state(_graph), _red_cost(_graph, _cost, _potential), kpeter@2581: _flow_result(0), _potential_result(0), kpeter@2581: _local_flow(false), _local_potential(false), kpeter@2575: _node_ref(graph), _edge_ref(graph) deba@2440: { kpeter@2575: // Copying _graph_ref to graph kpeter@2575: copyGraph(_graph, _graph_ref) kpeter@2575: .edgeMap(_capacity, capacity) kpeter@2575: .edgeMap(_cost, cost) kpeter@2575: .nodeRef(_node_ref) kpeter@2575: .edgeRef(_edge_ref) kpeter@2556: .run(); kpeter@2575: _supply[_node_ref[s]] = flow_value; kpeter@2575: _supply[_node_ref[t]] = -flow_value; kpeter@2575: _valid_supply = true; deba@2440: } deba@2440: kpeter@2581: /// Destructor. kpeter@2581: ~NetworkSimplex() { kpeter@2581: if (_local_flow) delete _flow_result; kpeter@2581: if (_local_potential) delete _potential_result; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Sets the flow map. kpeter@2581: /// kpeter@2581: /// Sets the flow map. kpeter@2581: /// kpeter@2581: /// \return \c (*this) kpeter@2581: NetworkSimplex& flowMap(FlowMap &map) { kpeter@2581: if (_local_flow) { kpeter@2581: delete _flow_result; kpeter@2581: _local_flow = false; kpeter@2581: } kpeter@2581: _flow_result = ↦ kpeter@2581: return *this; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Sets the potential map. kpeter@2581: /// kpeter@2581: /// Sets the potential map. kpeter@2581: /// kpeter@2581: /// \return \c (*this) kpeter@2581: NetworkSimplex& potentialMap(PotentialMap &map) { kpeter@2581: if (_local_potential) { kpeter@2581: delete _potential_result; kpeter@2581: _local_potential = false; kpeter@2581: } kpeter@2581: _potential_result = ↦ kpeter@2581: return *this; kpeter@2581: } kpeter@2581: kpeter@2581: /// \name Execution control kpeter@2581: /// The only way to execute the algorithm is to call the run() kpeter@2581: /// function. kpeter@2581: kpeter@2581: /// @{ kpeter@2581: kpeter@2556: /// \brief Runs the algorithm. kpeter@2556: /// kpeter@2556: /// Runs the algorithm. kpeter@2556: /// kpeter@2575: /// \param pivot_rule The pivot rule that is used during the kpeter@2575: /// algorithm. kpeter@2575: /// kpeter@2575: /// The available pivot rules: kpeter@2575: /// kpeter@2575: /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in kpeter@2575: /// a wraparound fashion in every iteration kpeter@2575: /// (\ref FirstEligiblePivotRule). kpeter@2575: /// kpeter@2575: /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in kpeter@2575: /// every iteration (\ref BestEligiblePivotRule). kpeter@2575: /// kpeter@2575: /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in kpeter@2575: /// every iteration in a wraparound fashion and the best eligible kpeter@2575: /// edge is selected from this block (\ref BlockSearchPivotRule). kpeter@2575: /// kpeter@2575: /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are kpeter@2575: /// examined in every iteration in a wraparound fashion and the best kpeter@2575: /// one is selected from them (\ref LimitedSearchPivotRule). kpeter@2575: /// kpeter@2575: /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is kpeter@2575: /// built from eligible edges and it is used for edge selection in kpeter@2575: /// the following minor iterations (\ref CandidateListPivotRule). kpeter@2575: /// kpeter@2575: /// - COMBINED_PIVOT This is a combined version of the two fastest kpeter@2575: /// pivot rules. kpeter@2575: /// For rather sparse graphs \ref LimitedSearchPivotRule kpeter@2575: /// "Limited Search" implementation is used, otherwise kpeter@2575: /// \ref BlockSearchPivotRule "Block Search" pivot rule is used. kpeter@2575: /// According to our benchmark tests this combined method is the kpeter@2575: /// most efficient. kpeter@2575: /// kpeter@2556: /// \return \c true if a feasible flow can be found. kpeter@2575: bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) { kpeter@2575: return init() && start(pivot_rule); kpeter@2556: } kpeter@2556: kpeter@2581: /// @} kpeter@2581: kpeter@2581: /// \name Query Functions kpeter@2581: /// The result of the algorithm can be obtained using these kpeter@2581: /// functions. kpeter@2581: /// \n run() must be called before using them. kpeter@2581: kpeter@2581: /// @{ kpeter@2581: kpeter@2575: /// \brief Returns a const reference to the edge map storing the kpeter@2575: /// found flow. deba@2440: /// kpeter@2575: /// Returns a const reference to the edge map storing the found flow. deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: const FlowMap& flowMap() const { kpeter@2581: return *_flow_result; deba@2440: } deba@2440: kpeter@2575: /// \brief Returns a const reference to the node map storing the kpeter@2575: /// found potentials (the dual solution). deba@2440: /// kpeter@2575: /// Returns a const reference to the node map storing the found kpeter@2575: /// potentials (the dual solution). deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: const PotentialMap& potentialMap() const { kpeter@2581: return *_potential_result; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Returns the flow on the edge. kpeter@2581: /// kpeter@2581: /// Returns the flow on the edge. kpeter@2581: /// kpeter@2581: /// \pre \ref run() must be called before using this function. kpeter@2581: Capacity flow(const typename Graph::Edge& edge) const { kpeter@2581: return (*_flow_result)[edge]; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Returns the potential of the node. kpeter@2581: /// kpeter@2581: /// Returns the potential of the node. kpeter@2581: /// kpeter@2581: /// \pre \ref run() must be called before using this function. kpeter@2581: Cost potential(const typename Graph::Node& node) const { kpeter@2581: return (*_potential_result)[node]; deba@2440: } deba@2440: deba@2440: /// \brief Returns the total cost of the found flow. deba@2440: /// deba@2440: /// Returns the total cost of the found flow. The complexity of the deba@2440: /// function is \f$ O(e) \f$. deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: Cost totalCost() const { deba@2440: Cost c = 0; kpeter@2575: for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) kpeter@2581: c += (*_flow_result)[e] * _cost[_edge_ref[e]]; deba@2440: return c; deba@2440: } deba@2440: kpeter@2581: /// @} kpeter@2581: kpeter@2575: private: deba@2440: deba@2440: /// \brief Extends the underlaying graph and initializes all the deba@2440: /// node and edge maps. deba@2440: bool init() { kpeter@2575: if (!_valid_supply) return false; deba@2440: kpeter@2581: // Initializing result maps kpeter@2581: if (!_flow_result) { kpeter@2581: _flow_result = new FlowMap(_graph_ref); kpeter@2581: _local_flow = true; kpeter@2581: } kpeter@2581: if (!_potential_result) { kpeter@2581: _potential_result = new PotentialMap(_graph_ref); kpeter@2581: _local_potential = true; kpeter@2581: } kpeter@2581: deba@2440: // Initializing state and flow maps kpeter@2575: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2575: _flow[e] = 0; kpeter@2575: _state[e] = STATE_LOWER; deba@2440: } deba@2440: deba@2440: // Adding an artificial root node to the graph kpeter@2575: _root = _graph.addNode(); kpeter@2575: _parent[_root] = INVALID; kpeter@2575: _pred_edge[_root] = INVALID; kpeter@2575: _depth[_root] = 0; kpeter@2575: _supply[_root] = 0; kpeter@2575: _potential[_root] = 0; deba@2440: deba@2440: // Adding artificial edges to the graph and initializing the node deba@2440: // maps of the spanning tree data structure kpeter@2575: Node last = _root; deba@2440: Edge e; deba@2440: Cost max_cost = std::numeric_limits::max() / 4; kpeter@2575: for (NodeIt u(_graph); u != INVALID; ++u) { kpeter@2575: if (u == _root) continue; kpeter@2575: _thread[last] = u; kpeter@2556: last = u; kpeter@2575: _parent[u] = _root; kpeter@2575: _depth[u] = 1; kpeter@2575: if (_supply[u] >= 0) { kpeter@2575: e = _graph.addEdge(u, _root); kpeter@2575: _flow[e] = _supply[u]; kpeter@2575: _forward[u] = true; kpeter@2575: _potential[u] = -max_cost; kpeter@2556: } else { kpeter@2575: e = _graph.addEdge(_root, u); kpeter@2575: _flow[e] = -_supply[u]; kpeter@2575: _forward[u] = false; kpeter@2575: _potential[u] = max_cost; kpeter@2556: } kpeter@2575: _cost[e] = max_cost; kpeter@2575: _capacity[e] = std::numeric_limits::max(); kpeter@2575: _state[e] = STATE_TREE; kpeter@2575: _pred_edge[u] = e; deba@2440: } kpeter@2575: _thread[last] = _root; deba@2440: kpeter@2575: return true; deba@2440: } deba@2440: kpeter@2575: /// Finds the join node. kpeter@2575: Node findJoinNode() { kpeter@2575: Node u = _graph.source(_in_edge); kpeter@2575: Node v = _graph.target(_in_edge); kpeter@2575: while (u != v) { kpeter@2575: if (_depth[u] == _depth[v]) { kpeter@2575: u = _parent[u]; kpeter@2575: v = _parent[v]; kpeter@2556: } kpeter@2575: else if (_depth[u] > _depth[v]) u = _parent[u]; kpeter@2575: else v = _parent[v]; deba@2440: } deba@2440: return u; deba@2440: } deba@2440: deba@2440: /// \brief Finds the leaving edge of the cycle. Returns \c true if deba@2440: /// the leaving edge is not the same as the entering edge. deba@2440: bool findLeavingEdge() { deba@2440: // Initializing first and second nodes according to the direction deba@2440: // of the cycle kpeter@2575: if (_state[_in_edge] == STATE_LOWER) { kpeter@2575: first = _graph.source(_in_edge); kpeter@2575: second = _graph.target(_in_edge); deba@2440: } else { kpeter@2575: first = _graph.target(_in_edge); kpeter@2575: second = _graph.source(_in_edge); deba@2440: } kpeter@2575: delta = _capacity[_in_edge]; deba@2440: bool result = false; deba@2440: Capacity d; deba@2440: Edge e; deba@2440: deba@2440: // Searching the cycle along the path form the first node to the deba@2440: // root node kpeter@2575: for (Node u = first; u != join; u = _parent[u]) { kpeter@2575: e = _pred_edge[u]; kpeter@2575: d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e]; kpeter@2556: if (d < delta) { kpeter@2556: delta = d; kpeter@2556: u_out = u; kpeter@2556: u_in = first; kpeter@2556: v_in = second; kpeter@2556: result = true; kpeter@2556: } deba@2440: } deba@2440: // Searching the cycle along the path form the second node to the deba@2440: // root node kpeter@2575: for (Node u = second; u != join; u = _parent[u]) { kpeter@2575: e = _pred_edge[u]; kpeter@2575: d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e]; kpeter@2556: if (d <= delta) { kpeter@2556: delta = d; kpeter@2556: u_out = u; kpeter@2556: u_in = second; kpeter@2556: v_in = first; kpeter@2556: result = true; kpeter@2556: } deba@2440: } deba@2440: return result; deba@2440: } deba@2440: kpeter@2575: /// Changes \c flow and \c state edge maps. deba@2440: void changeFlows(bool change) { deba@2440: // Augmenting along the cycle deba@2440: if (delta > 0) { kpeter@2575: Capacity val = _state[_in_edge] * delta; kpeter@2575: _flow[_in_edge] += val; kpeter@2575: for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) { kpeter@2575: _flow[_pred_edge[u]] += _forward[u] ? -val : val; kpeter@2556: } kpeter@2575: for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) { kpeter@2575: _flow[_pred_edge[u]] += _forward[u] ? val : -val; kpeter@2556: } deba@2440: } deba@2440: // Updating the state of the entering and leaving edges deba@2440: if (change) { kpeter@2575: _state[_in_edge] = STATE_TREE; kpeter@2575: _state[_pred_edge[u_out]] = kpeter@2575: (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; deba@2440: } else { kpeter@2575: _state[_in_edge] = -_state[_in_edge]; deba@2440: } deba@2440: } deba@2440: kpeter@2575: /// Updates \c thread and \c parent node maps. deba@2440: void updateThreadParent() { deba@2440: Node u; kpeter@2575: v_out = _parent[u_out]; deba@2440: deba@2440: // Handling the case when join and v_out coincide deba@2440: bool par_first = false; deba@2440: if (join == v_out) { kpeter@2575: for (u = join; u != u_in && u != v_in; u = _thread[u]) ; kpeter@2556: if (u == v_in) { kpeter@2556: par_first = true; kpeter@2575: while (_thread[u] != u_out) u = _thread[u]; kpeter@2556: first = u; kpeter@2556: } deba@2440: } deba@2440: deba@2440: // Finding the last successor of u_in (u) and the node after it deba@2440: // (right) according to the thread index kpeter@2575: for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ; kpeter@2575: right = _thread[u]; kpeter@2575: if (_thread[v_in] == u_out) { kpeter@2575: for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ; kpeter@2575: if (last == u_out) last = _thread[last]; deba@2440: } kpeter@2575: else last = _thread[v_in]; deba@2440: deba@2440: // Updating stem nodes kpeter@2575: _thread[v_in] = stem = u_in; deba@2440: par_stem = v_in; deba@2440: while (stem != u_out) { kpeter@2575: _thread[u] = new_stem = _parent[stem]; deba@2440: kpeter@2556: // Finding the node just before the stem node (u) according to kpeter@2556: // the original thread index kpeter@2575: for (u = new_stem; _thread[u] != stem; u = _thread[u]) ; kpeter@2575: _thread[u] = right; deba@2440: kpeter@2556: // Changing the parent node of stem and shifting stem and kpeter@2556: // par_stem nodes kpeter@2575: _parent[stem] = par_stem; kpeter@2556: par_stem = stem; kpeter@2556: stem = new_stem; deba@2440: kpeter@2556: // Finding the last successor of stem (u) and the node after it kpeter@2556: // (right) according to the thread index kpeter@2575: for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ; kpeter@2575: right = _thread[u]; deba@2440: } kpeter@2575: _parent[u_out] = par_stem; kpeter@2575: _thread[u] = last; deba@2440: deba@2440: if (join == v_out && par_first) { kpeter@2575: if (first != v_in) _thread[first] = right; deba@2440: } else { kpeter@2575: for (u = v_out; _thread[u] != u_out; u = _thread[u]) ; kpeter@2575: _thread[u] = right; deba@2440: } deba@2440: } deba@2440: kpeter@2575: /// Updates \c pred_edge and \c forward node maps. deba@2440: void updatePredEdge() { deba@2440: Node u = u_out, v; deba@2440: while (u != u_in) { kpeter@2575: v = _parent[u]; kpeter@2575: _pred_edge[u] = _pred_edge[v]; kpeter@2575: _forward[u] = !_forward[v]; kpeter@2556: u = v; deba@2440: } kpeter@2575: _pred_edge[u_in] = _in_edge; kpeter@2575: _forward[u_in] = (u_in == _graph.source(_in_edge)); deba@2440: } deba@2440: kpeter@2575: /// Updates \c depth and \c potential node maps. deba@2440: void updateDepthPotential() { kpeter@2575: _depth[u_in] = _depth[v_in] + 1; kpeter@2575: _potential[u_in] = _forward[u_in] ? kpeter@2575: _potential[v_in] - _cost[_pred_edge[u_in]] : kpeter@2575: _potential[v_in] + _cost[_pred_edge[u_in]]; deba@2440: kpeter@2575: Node u = _thread[u_in], v; deba@2440: while (true) { kpeter@2575: v = _parent[u]; kpeter@2556: if (v == INVALID) break; kpeter@2575: _depth[u] = _depth[v] + 1; kpeter@2575: _potential[u] = _forward[u] ? kpeter@2575: _potential[v] - _cost[_pred_edge[u]] : kpeter@2575: _potential[v] + _cost[_pred_edge[u]]; kpeter@2575: if (_depth[u] <= _depth[v_in]) break; kpeter@2575: u = _thread[u]; deba@2440: } deba@2440: } deba@2440: kpeter@2575: /// Executes the algorithm. kpeter@2575: bool start(PivotRuleEnum pivot_rule) { kpeter@2575: switch (pivot_rule) { kpeter@2575: case FIRST_ELIGIBLE_PIVOT: kpeter@2575: return start(); kpeter@2575: case BEST_ELIGIBLE_PIVOT: kpeter@2575: return start(); kpeter@2575: case BLOCK_SEARCH_PIVOT: kpeter@2575: return start(); kpeter@2575: case LIMITED_SEARCH_PIVOT: kpeter@2575: return start(); kpeter@2575: case CANDIDATE_LIST_PIVOT: kpeter@2575: return start(); kpeter@2575: case COMBINED_PIVOT: kpeter@2575: if ( countEdges(_graph) / countNodes(_graph) <= kpeter@2575: COMBINED_PIVOT_MAX_DEG ) kpeter@2575: return start(); kpeter@2575: else kpeter@2575: return start(); kpeter@2575: } kpeter@2575: return false; kpeter@2575: } kpeter@2575: kpeter@2575: template deba@2440: bool start() { kpeter@2575: PivotRuleImplementation pivot(*this); kpeter@2575: kpeter@2575: // Executing the network simplex algorithm kpeter@2575: while (pivot.findEnteringEdge()) { kpeter@2556: join = findJoinNode(); kpeter@2556: bool change = findLeavingEdge(); kpeter@2556: changeFlows(change); kpeter@2556: if (change) { kpeter@2556: updateThreadParent(); kpeter@2556: updatePredEdge(); kpeter@2556: updateDepthPotential(); kpeter@2556: } deba@2440: } deba@2440: kpeter@2575: // Checking if the flow amount equals zero on all the artificial kpeter@2575: // edges kpeter@2575: for (InEdgeIt e(_graph, _root); e != INVALID; ++e) kpeter@2575: if (_flow[e] > 0) return false; kpeter@2575: for (OutEdgeIt e(_graph, _root); e != INVALID; ++e) kpeter@2575: if (_flow[e] > 0) return false; deba@2440: kpeter@2575: // Copying flow values to _flow_result kpeter@2575: if (_lower) { kpeter@2575: for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) kpeter@2581: (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]]; deba@2440: } else { kpeter@2575: for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) kpeter@2581: (*_flow_result)[e] = _flow[_edge_ref[e]]; deba@2440: } kpeter@2575: // Copying potential values to _potential_result kpeter@2575: for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) kpeter@2581: (*_potential_result)[n] = _potential[_node_ref[n]]; deba@2440: deba@2440: return true; deba@2440: } deba@2440: deba@2440: }; //class NetworkSimplex deba@2440: deba@2440: ///@} deba@2440: deba@2440: } //namespace lemon deba@2440: deba@2440: #endif //LEMON_NETWORK_SIMPLEX_H