deba@2440: /* -*- C++ -*- deba@2440: * deba@2440: * This file is a part of LEMON, a generic C++ optimization library deba@2440: * deba@2440: * Copyright (C) 2003-2007 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 deba@2440: /// \brief The network simplex algorithm for finding a minimum cost deba@2440: /// flow. deba@2440: deba@2440: #include kpeter@2509: #include kpeter@2509: #include deba@2440: #include deba@2440: deba@2440: /// \brief The pivot rule used in the algorithm. deba@2440: //#define FIRST_ELIGIBLE_PIVOT deba@2440: //#define BEST_ELIGIBLE_PIVOT deba@2444: #define EDGE_BLOCK_PIVOT deba@2440: //#define CANDIDATE_LIST_PIVOT deba@2440: //#define SORTED_LIST_PIVOT deba@2440: deba@2444: //#define _DEBUG_ITER_ deba@2444: deba@2444: deba@2440: /// \brief State constant for edges at their lower bounds. deba@2440: #define LOWER 1 deba@2440: /// \brief State constant for edges in the spanning tree. deba@2440: #define TREE 0 deba@2440: /// \brief State constant for edges at their upper bounds. deba@2440: #define UPPER -1 deba@2440: deba@2440: #ifdef EDGE_BLOCK_PIVOT kpeter@2471: #include deba@2444: /// \brief Lower bound for the size of blocks. deba@2444: #define MIN_BLOCK_SIZE 10 deba@2440: #endif deba@2440: deba@2440: #ifdef CANDIDATE_LIST_PIVOT kpeter@2471: #include kpeter@2471: #define LIST_LENGTH_DIV 50 kpeter@2471: #define MINOR_LIMIT_DIV 200 deba@2440: #endif deba@2440: deba@2440: #ifdef SORTED_LIST_PIVOT kpeter@2471: #include deba@2440: #include kpeter@2471: #define LIST_LENGTH_DIV 100 kpeter@2471: #define LOWER_DIV 4 deba@2440: #endif 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: /// deba@2440: /// \ref lemon::NetworkSimplex "NetworkSimplex" implements the deba@2440: /// network simplex algorithm for finding a minimum cost flow. deba@2440: /// deba@2440: /// \param Graph The directed graph type the algorithm runs on. deba@2440: /// \param LowerMap The type of the lower bound map. deba@2440: /// \param CapacityMap The type of the capacity (upper bound) map. deba@2440: /// \param CostMap The type of the cost (length) map. deba@2440: /// \param SupplyMap The type of the supply map. deba@2440: /// deba@2440: /// \warning deba@2440: /// - Edge capacities and costs should be nonnegative integers. deba@2440: /// However \c CostMap::Value should be signed type. kpeter@2509: /// - Supply values should be signed integers. deba@2440: /// - \c LowerMap::Value must be convertible to deba@2440: /// \c CapacityMap::Value and \c CapacityMap::Value must be deba@2440: /// convertible to \c SupplyMap::Value. deba@2440: /// deba@2440: /// \author Peter Kovacs deba@2440: kpeter@2533: template < typename Graph, kpeter@2533: typename LowerMap = typename Graph::template EdgeMap, kpeter@2533: typename CapacityMap = LowerMap, kpeter@2533: typename CostMap = typename Graph::template EdgeMap, kpeter@2533: typename SupplyMap = typename Graph::template NodeMap kpeter@2533: > deba@2440: class NetworkSimplex deba@2440: { deba@2440: typedef typename LowerMap::Value Lower; 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; deba@2440: typedef typename SGraph::Node Node; deba@2440: typedef typename SGraph::NodeIt NodeIt; deba@2440: typedef typename SGraph::Edge Edge; deba@2440: typedef typename SGraph::EdgeIt EdgeIt; deba@2440: typedef typename SGraph::InEdgeIt InEdgeIt; deba@2440: typedef typename SGraph::OutEdgeIt OutEdgeIt; deba@2440: deba@2440: typedef typename SGraph::template EdgeMap SLowerMap; 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: deba@2440: /// \brief The type of the flow map. deba@2440: typedef typename Graph::template EdgeMap FlowMap; deba@2440: /// \brief The type of the potential map. deba@2440: typedef typename Graph::template NodeMap PotentialMap; deba@2440: deba@2440: protected: deba@2440: deba@2440: /// \brief Map adaptor class for handling reduced edge costs. deba@2440: class ReducedCostMap : public MapBase deba@2440: { deba@2440: private: deba@2440: deba@2440: const SGraph &gr; deba@2440: const SCostMap &cost_map; deba@2440: const SPotentialMap &pot_map; deba@2440: deba@2440: public: deba@2440: deba@2440: ReducedCostMap( const SGraph &_gr, deba@2440: const SCostMap &_cm, deba@2440: const SPotentialMap &_pm ) : deba@2440: gr(_gr), cost_map(_cm), pot_map(_pm) {} deba@2440: kpeter@2509: Cost operator[](const Edge &e) const { deba@2440: return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)]; deba@2440: } deba@2440: deba@2440: }; //class ReducedCostMap deba@2440: deba@2440: protected: deba@2440: deba@2440: /// \brief The directed graph the algorithm runs on. deba@2440: SGraph graph; deba@2440: /// \brief The original graph. deba@2440: const Graph &graph_ref; deba@2440: /// \brief The original lower bound map. deba@2440: const LowerMap *lower; deba@2440: /// \brief The capacity map. deba@2440: SCapacityMap capacity; deba@2440: /// \brief The cost map. deba@2440: SCostMap cost; deba@2440: /// \brief The supply map. deba@2440: SSupplyMap supply; deba@2440: /// \brief The reduced cost map. deba@2440: ReducedCostMap red_cost; deba@2440: /// \brief The sum of supply values equals zero. deba@2440: bool valid_supply; deba@2440: deba@2440: /// \brief The edge map of the current flow. deba@2440: SCapacityMap flow; deba@2440: /// \brief The edge map of the found flow on the original graph. deba@2440: FlowMap flow_result; deba@2440: /// \brief The potential node map. deba@2440: SPotentialMap potential; deba@2440: /// \brief The potential node map on the original graph. deba@2440: PotentialMap potential_result; deba@2440: deba@2440: /// \brief Node reference for the original graph. deba@2440: NodeRefMap node_ref; deba@2440: /// \brief Edge reference for the original graph. deba@2440: EdgeRefMap edge_ref; deba@2440: deba@2440: /// \brief The depth node map of the spanning tree structure. deba@2440: IntNodeMap depth; deba@2440: /// \brief The parent node map of the spanning tree structure. deba@2440: NodeNodeMap parent; deba@2440: /// \brief The pred_edge node map of the spanning tree structure. deba@2440: EdgeNodeMap pred_edge; deba@2440: /// \brief The thread node map of the spanning tree structure. deba@2440: NodeNodeMap thread; deba@2440: /// \brief The forward node map of the spanning tree structure. deba@2440: BoolNodeMap forward; deba@2440: /// \brief The state edge map. deba@2440: IntEdgeMap state; deba@2440: deba@2440: deba@2440: #ifdef EDGE_BLOCK_PIVOT deba@2440: /// \brief The size of blocks for the "Edge Block" pivot rule. deba@2440: int block_size; deba@2440: #endif deba@2440: #ifdef CANDIDATE_LIST_PIVOT deba@2440: /// \brief The list of candidate edges for the "Candidate List" deba@2440: /// pivot rule. kpeter@2471: std::vector candidates; kpeter@2471: /// \brief The maximum length of the edge list for the kpeter@2471: /// "Candidate List" pivot rule. kpeter@2471: int list_length; kpeter@2471: /// \brief The maximum number of minor iterations between two major kpeter@2471: /// itarations. kpeter@2471: int minor_limit; deba@2440: /// \brief The number of minor iterations. deba@2440: int minor_count; deba@2440: #endif deba@2440: #ifdef SORTED_LIST_PIVOT deba@2440: /// \brief The list of candidate edges for the "Sorted List" deba@2440: /// pivot rule. kpeter@2471: std::vector candidates; kpeter@2471: /// \brief The maximum length of the edge list for the kpeter@2471: /// "Sorted List" pivot rule. kpeter@2471: int list_length; kpeter@2471: int list_index; deba@2440: #endif deba@2440: deba@2440: // Root node of the starting spanning tree. deba@2440: Node root; deba@2440: // The entering edge of the current pivot iteration. deba@2440: Edge in_edge; deba@2440: // Temporary nodes used in the current pivot iteration. deba@2440: Node join, u_in, v_in, u_out, v_out; deba@2440: Node right, first, second, last; deba@2440: Node stem, par_stem, new_stem; deba@2440: // The maximum augment amount along the cycle in the current pivot deba@2440: // iteration. deba@2440: Capacity delta; deba@2440: deba@2440: public : deba@2440: deba@2440: /// \brief General constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// General constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _lower The lower bounds of the edges. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _supply The supply values of the nodes (signed). deba@2440: NetworkSimplex( const Graph &_graph, deba@2440: const LowerMap &_lower, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: const SupplyMap &_supply ) : deba@2440: graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), deba@2440: supply(graph), flow(graph), flow_result(_graph), potential(graph), deba@2440: potential_result(_graph), depth(graph), parent(graph), deba@2440: pred_edge(graph), thread(graph), forward(graph), state(graph), deba@2440: node_ref(graph_ref), edge_ref(graph_ref), deba@2440: red_cost(graph, cost, potential) deba@2440: { deba@2440: // Checking the sum of supply values deba@2440: Supply sum = 0; deba@2440: for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) deba@2440: sum += _supply[n]; deba@2440: if (!(valid_supply = sum == 0)) return; deba@2440: deba@2440: // Copying graph_ref to graph deba@2457: graph.reserveNode(countNodes(graph_ref) + 1); deba@2457: graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref)); deba@2440: copyGraph(graph, graph_ref) deba@2440: .edgeMap(cost, _cost) deba@2440: .nodeRef(node_ref) deba@2440: .edgeRef(edge_ref) deba@2440: .run(); deba@2440: deba@2440: // Removing nonzero lower bounds deba@2440: for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { deba@2440: capacity[edge_ref[e]] = _capacity[e] - _lower[e]; deba@2440: } deba@2440: for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { deba@2440: Supply s = _supply[n]; deba@2440: for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) deba@2440: s += _lower[e]; deba@2440: for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) deba@2440: s -= _lower[e]; deba@2440: supply[node_ref[n]] = s; deba@2440: } deba@2440: } deba@2440: deba@2440: /// \brief General constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// General constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _supply The supply values of the nodes (signed). deba@2440: NetworkSimplex( const Graph &_graph, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: const SupplyMap &_supply ) : deba@2440: graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), deba@2440: supply(graph), flow(graph), flow_result(_graph), potential(graph), deba@2440: potential_result(_graph), depth(graph), parent(graph), deba@2440: pred_edge(graph), thread(graph), forward(graph), state(graph), deba@2440: node_ref(graph_ref), edge_ref(graph_ref), deba@2440: red_cost(graph, cost, potential) deba@2440: { deba@2440: // Checking the sum of supply values deba@2440: Supply sum = 0; deba@2440: for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) deba@2440: sum += _supply[n]; deba@2440: if (!(valid_supply = sum == 0)) return; deba@2440: deba@2440: // Copying graph_ref to graph deba@2440: copyGraph(graph, graph_ref) deba@2440: .edgeMap(capacity, _capacity) deba@2440: .edgeMap(cost, _cost) deba@2440: .nodeMap(supply, _supply) deba@2440: .nodeRef(node_ref) deba@2440: .edgeRef(edge_ref) deba@2440: .run(); deba@2440: } deba@2440: deba@2440: /// \brief Simple constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// Simple constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _lower The lower bounds of the edges. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _s The source node. deba@2440: /// \param _t The target node. deba@2440: /// \param _flow_value The required amount of flow from node \c _s deba@2440: /// to node \c _t (i.e. the supply of \c _s and the demand of deba@2440: /// \c _t). deba@2440: NetworkSimplex( const Graph &_graph, deba@2440: const LowerMap &_lower, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: typename Graph::Node _s, deba@2440: typename Graph::Node _t, deba@2440: typename SupplyMap::Value _flow_value ) : deba@2440: graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), deba@2440: supply(graph), flow(graph), flow_result(_graph), potential(graph), deba@2440: potential_result(_graph), depth(graph), parent(graph), deba@2440: pred_edge(graph), thread(graph), forward(graph), state(graph), deba@2440: node_ref(graph_ref), edge_ref(graph_ref), deba@2440: red_cost(graph, cost, potential) deba@2440: { deba@2440: // Copying graph_ref to graph deba@2440: copyGraph(graph, graph_ref) deba@2440: .edgeMap(cost, _cost) deba@2440: .nodeRef(node_ref) deba@2440: .edgeRef(edge_ref) deba@2440: .run(); deba@2440: deba@2440: // Removing nonzero lower bounds deba@2440: for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { deba@2440: capacity[edge_ref[e]] = _capacity[e] - _lower[e]; deba@2440: } deba@2440: for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { deba@2440: Supply s = 0; deba@2440: if (n == _s) s = _flow_value; deba@2440: if (n == _t) s = -_flow_value; deba@2440: for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) deba@2440: s += _lower[e]; deba@2440: for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) deba@2440: s -= _lower[e]; deba@2440: supply[node_ref[n]] = s; deba@2440: } deba@2440: valid_supply = true; deba@2440: } deba@2440: deba@2440: /// \brief Simple constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// Simple constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _s The source node. deba@2440: /// \param _t The target node. deba@2440: /// \param _flow_value The required amount of flow from node \c _s deba@2440: /// to node \c _t (i.e. the supply of \c _s and the demand of deba@2440: /// \c _t). deba@2440: NetworkSimplex( const Graph &_graph, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: typename Graph::Node _s, deba@2440: typename Graph::Node _t, deba@2440: typename SupplyMap::Value _flow_value ) : deba@2440: graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), deba@2440: supply(graph, 0), flow(graph), flow_result(_graph), potential(graph), deba@2440: potential_result(_graph), depth(graph), parent(graph), deba@2440: pred_edge(graph), thread(graph), forward(graph), state(graph), deba@2440: node_ref(graph_ref), edge_ref(graph_ref), deba@2440: red_cost(graph, cost, potential) deba@2440: { deba@2440: // Copying graph_ref to graph deba@2440: copyGraph(graph, graph_ref) deba@2440: .edgeMap(capacity, _capacity) deba@2440: .edgeMap(cost, _cost) deba@2440: .nodeRef(node_ref) deba@2440: .edgeRef(edge_ref) deba@2440: .run(); deba@2440: supply[node_ref[_s]] = _flow_value; deba@2440: supply[node_ref[_t]] = -_flow_value; deba@2440: valid_supply = true; deba@2440: } deba@2440: deba@2440: /// \brief Returns a const reference to the flow map. deba@2440: /// deba@2440: /// Returns a const reference to the flow map. deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: const FlowMap& flowMap() const { deba@2440: return flow_result; deba@2440: } deba@2440: deba@2440: /// \brief Returns a const reference to the potential map (the dual deba@2440: /// solution). deba@2440: /// deba@2440: /// Returns a const reference to the potential map (the dual deba@2440: /// solution). deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: const PotentialMap& potentialMap() const { deba@2440: return potential_result; 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; deba@2440: for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) deba@2440: c += flow_result[e] * cost[edge_ref[e]]; deba@2440: return c; deba@2440: } deba@2440: deba@2440: /// \brief Runs the algorithm. deba@2440: /// deba@2440: /// Runs the algorithm. deba@2440: /// deba@2440: /// \return \c true if a feasible flow can be found. deba@2440: bool run() { deba@2440: return init() && start(); deba@2440: } deba@2440: deba@2440: protected: deba@2440: deba@2440: /// \brief Extends the underlaying graph and initializes all the deba@2440: /// node and edge maps. deba@2440: bool init() { deba@2440: if (!valid_supply) return false; deba@2440: deba@2440: // Initializing state and flow maps deba@2440: for (EdgeIt e(graph); e != INVALID; ++e) { deba@2440: flow[e] = 0; deba@2440: state[e] = LOWER; deba@2440: } deba@2440: deba@2440: // Adding an artificial root node to the graph deba@2440: root = graph.addNode(); deba@2440: parent[root] = INVALID; deba@2440: pred_edge[root] = INVALID; deba@2457: depth[root] = 0; deba@2457: supply[root] = 0; deba@2457: 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 deba@2440: Supply sum = 0; deba@2440: Node last = root; deba@2440: Edge e; deba@2440: Cost max_cost = std::numeric_limits::max() / 4; deba@2440: for (NodeIt u(graph); u != INVALID; ++u) { deba@2440: if (u == root) continue; deba@2440: thread[last] = u; deba@2440: last = u; deba@2440: parent[u] = root; deba@2440: depth[u] = 1; deba@2440: sum += supply[u]; deba@2440: if (supply[u] >= 0) { deba@2440: e = graph.addEdge(u, root); deba@2440: flow[e] = supply[u]; deba@2440: forward[u] = true; deba@2440: potential[u] = max_cost; deba@2440: } else { deba@2440: e = graph.addEdge(root, u); deba@2440: flow[e] = -supply[u]; deba@2440: forward[u] = false; deba@2440: potential[u] = -max_cost; deba@2440: } deba@2440: cost[e] = max_cost; deba@2440: capacity[e] = std::numeric_limits::max(); deba@2440: state[e] = TREE; deba@2440: pred_edge[u] = e; deba@2440: } deba@2440: thread[last] = root; deba@2440: deba@2440: #ifdef EDGE_BLOCK_PIVOT deba@2440: // Initializing block_size for the edge block pivot rule deba@2440: int edge_num = countEdges(graph); kpeter@2471: block_size = 2 * int(sqrt(countEdges(graph))); kpeter@2471: if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE; kpeter@2471: // block_size = edge_num >= BLOCK_NUM * MIN_BLOCK_SIZE ? kpeter@2471: // edge_num / BLOCK_NUM : MIN_BLOCK_SIZE; deba@2440: #endif deba@2440: #ifdef CANDIDATE_LIST_PIVOT kpeter@2471: int edge_num = countEdges(graph); deba@2440: minor_count = 0; kpeter@2471: list_length = edge_num / LIST_LENGTH_DIV; kpeter@2471: minor_limit = edge_num / MINOR_LIMIT_DIV; kpeter@2471: #endif kpeter@2471: #ifdef SORTED_LIST_PIVOT kpeter@2471: int edge_num = countEdges(graph); kpeter@2471: list_index = 0; kpeter@2471: list_length = edge_num / LIST_LENGTH_DIV; deba@2440: #endif deba@2440: deba@2440: return sum == 0; deba@2440: } deba@2440: deba@2440: #ifdef FIRST_ELIGIBLE_PIVOT deba@2440: /// \brief Finds entering edge according to the "First Eligible" deba@2440: /// pivot rule. deba@2440: bool findEnteringEdge(EdgeIt &next_edge) { deba@2440: for (EdgeIt e = next_edge; e != INVALID; ++e) { deba@2440: if (state[e] * red_cost[e] < 0) { deba@2440: in_edge = e; deba@2440: next_edge = ++e; deba@2440: return true; deba@2440: } deba@2440: } deba@2440: for (EdgeIt e(graph); e != next_edge; ++e) { deba@2440: if (state[e] * red_cost[e] < 0) { deba@2440: in_edge = e; deba@2440: next_edge = ++e; deba@2440: return true; deba@2440: } deba@2440: } deba@2440: return false; deba@2440: } deba@2440: #endif deba@2440: deba@2440: #ifdef BEST_ELIGIBLE_PIVOT deba@2440: /// \brief Finds entering edge according to the "Best Eligible" deba@2440: /// pivot rule. deba@2440: bool findEnteringEdge() { deba@2440: Cost min = 0; deba@2440: for (EdgeIt e(graph); e != INVALID; ++e) { deba@2440: if (state[e] * red_cost[e] < min) { deba@2440: min = state[e] * red_cost[e]; deba@2440: in_edge = e; deba@2440: } deba@2440: } deba@2440: return min < 0; deba@2440: } deba@2440: #endif deba@2440: deba@2440: #ifdef EDGE_BLOCK_PIVOT deba@2440: /// \brief Finds entering edge according to the "Edge Block" deba@2440: /// pivot rule. deba@2440: bool findEnteringEdge(EdgeIt &next_edge) { deba@2444: // Performing edge block selection deba@2444: Cost curr, min = 0; deba@2444: EdgeIt min_edge(graph); deba@2444: int cnt = 0; deba@2444: for (EdgeIt e = next_edge; e != INVALID; ++e) { deba@2444: if ((curr = state[e] * red_cost[e]) < min) { deba@2444: min = curr; deba@2444: min_edge = e; deba@2440: } deba@2444: if (++cnt == block_size) { deba@2444: if (min < 0) break; deba@2444: cnt = 0; deba@2444: } deba@2444: } deba@2444: if (!(min < 0)) { deba@2440: for (EdgeIt e(graph); e != next_edge; ++e) { deba@2440: if ((curr = state[e] * red_cost[e]) < min) { deba@2440: min = curr; deba@2440: min_edge = e; deba@2440: } deba@2440: if (++cnt == block_size) { deba@2440: if (min < 0) break; deba@2440: cnt = 0; deba@2440: } deba@2440: } deba@2440: } deba@2444: in_edge = min_edge; deba@2444: if ((next_edge = ++min_edge) == INVALID) deba@2444: next_edge = EdgeIt(graph); deba@2444: return min < 0; deba@2440: } deba@2440: #endif deba@2440: deba@2440: #ifdef CANDIDATE_LIST_PIVOT deba@2440: /// \brief Finds entering edge according to the "Candidate List" deba@2440: /// pivot rule. deba@2440: bool findEnteringEdge() { kpeter@2471: typedef typename std::vector::iterator ListIt; deba@2440: kpeter@2471: if (minor_count >= minor_limit || candidates.size() == 0) { deba@2440: // Major iteration kpeter@2471: candidates.clear(); deba@2440: for (EdgeIt e(graph); e != INVALID; ++e) { deba@2440: if (state[e] * red_cost[e] < 0) { deba@2440: candidates.push_back(e); kpeter@2471: if (candidates.size() == list_length) break; deba@2440: } deba@2440: } deba@2440: if (candidates.size() == 0) return false; deba@2440: } deba@2440: deba@2440: // Minor iteration deba@2440: ++minor_count; deba@2440: Cost min = 0; kpeter@2471: Edge e; kpeter@2471: for (int i = 0; i < candidates.size(); ++i) { kpeter@2471: e = candidates[i]; kpeter@2471: if (state[e] * red_cost[e] < min) { kpeter@2471: min = state[e] * red_cost[e]; kpeter@2471: in_edge = e; deba@2440: } deba@2440: } deba@2440: return true; deba@2440: } deba@2440: #endif deba@2440: deba@2440: #ifdef SORTED_LIST_PIVOT deba@2440: /// \brief Functor class to compare edges during sort of the deba@2440: /// candidate list. deba@2440: class SortFunc deba@2440: { deba@2440: private: deba@2440: const IntEdgeMap &st; deba@2440: const ReducedCostMap &rc; deba@2440: public: deba@2440: SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) : deba@2440: st(_st), rc(_rc) {} deba@2440: bool operator()(const Edge &e1, const Edge &e2) { deba@2440: return st[e1] * rc[e1] < st[e2] * rc[e2]; deba@2440: } deba@2440: }; deba@2440: deba@2440: /// \brief Finds entering edge according to the "Sorted List" deba@2440: /// pivot rule. deba@2440: bool findEnteringEdge() { deba@2440: static SortFunc sort_func(state, red_cost); deba@2440: deba@2440: // Minor iteration kpeter@2471: while (list_index < candidates.size()) { kpeter@2471: in_edge = candidates[list_index++]; deba@2440: if (state[in_edge] * red_cost[in_edge] < 0) return true; deba@2440: } deba@2440: deba@2440: // Major iteration kpeter@2471: candidates.clear(); deba@2440: Cost curr, min = 0; deba@2440: for (EdgeIt e(graph); e != INVALID; ++e) { deba@2440: if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) { deba@2440: candidates.push_back(e); deba@2440: if (curr < min) min = curr; kpeter@2471: if (candidates.size() == list_length) break; deba@2440: } deba@2440: } deba@2440: if (candidates.size() == 0) return false; deba@2440: sort(candidates.begin(), candidates.end(), sort_func); kpeter@2471: in_edge = candidates[0]; kpeter@2471: list_index = 1; deba@2440: return true; deba@2440: } deba@2440: #endif deba@2440: deba@2440: /// \brief Finds the join node. deba@2440: Node findJoinNode() { deba@2440: // Finding the join node deba@2440: Node u = graph.source(in_edge); deba@2440: Node v = graph.target(in_edge); deba@2440: while (u != v) { deba@2440: if (depth[u] == depth[v]) { deba@2440: u = parent[u]; deba@2440: v = parent[v]; deba@2440: } deba@2440: else if (depth[u] > depth[v]) u = parent[u]; deba@2440: 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 deba@2440: if (state[in_edge] == LOWER) { deba@2440: first = graph.source(in_edge); deba@2440: second = graph.target(in_edge); deba@2440: } else { deba@2440: first = graph.target(in_edge); deba@2440: second = graph.source(in_edge); deba@2440: } deba@2440: 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 deba@2440: for (Node u = first; u != join; u = parent[u]) { deba@2440: e = pred_edge[u]; deba@2440: d = forward[u] ? flow[e] : capacity[e] - flow[e]; deba@2440: if (d < delta) { deba@2440: delta = d; deba@2440: u_out = u; deba@2440: u_in = first; deba@2440: v_in = second; deba@2440: result = true; deba@2440: } deba@2440: } deba@2440: // Searching the cycle along the path form the second node to the deba@2440: // root node deba@2440: for (Node u = second; u != join; u = parent[u]) { deba@2440: e = pred_edge[u]; deba@2440: d = forward[u] ? capacity[e] - flow[e] : flow[e]; deba@2440: if (d <= delta) { deba@2440: delta = d; deba@2440: u_out = u; deba@2440: u_in = second; deba@2440: v_in = first; deba@2440: result = true; deba@2440: } deba@2440: } deba@2440: return result; deba@2440: } deba@2440: deba@2440: /// \brief Changes flow and state edge maps. deba@2440: void changeFlows(bool change) { deba@2440: // Augmenting along the cycle deba@2440: if (delta > 0) { deba@2440: Capacity val = state[in_edge] * delta; deba@2440: flow[in_edge] += val; deba@2440: for (Node u = graph.source(in_edge); u != join; u = parent[u]) { deba@2440: flow[pred_edge[u]] += forward[u] ? -val : val; deba@2440: } deba@2440: for (Node u = graph.target(in_edge); u != join; u = parent[u]) { deba@2440: flow[pred_edge[u]] += forward[u] ? val : -val; deba@2440: } deba@2440: } deba@2440: // Updating the state of the entering and leaving edges deba@2440: if (change) { deba@2440: state[in_edge] = TREE; deba@2440: state[pred_edge[u_out]] = deba@2440: (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER; deba@2440: } else { deba@2440: state[in_edge] = -state[in_edge]; deba@2440: } deba@2440: } deba@2440: deba@2440: /// \brief Updates thread and parent node maps. deba@2440: void updateThreadParent() { deba@2440: Node u; deba@2440: 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) { deba@2440: for (u = join; u != u_in && u != v_in; u = thread[u]) ; deba@2440: if (u == v_in) { deba@2440: par_first = true; deba@2440: while (thread[u] != u_out) u = thread[u]; deba@2440: first = u; deba@2440: } 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 deba@2440: for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ; deba@2440: right = thread[u]; deba@2440: if (thread[v_in] == u_out) { deba@2440: for (last = u; depth[last] > depth[u_out]; last = thread[last]) ; deba@2440: if (last == u_out) last = thread[last]; deba@2440: } deba@2440: else last = thread[v_in]; deba@2440: deba@2440: // Updating stem nodes deba@2440: thread[v_in] = stem = u_in; deba@2440: par_stem = v_in; deba@2440: while (stem != u_out) { deba@2440: thread[u] = new_stem = parent[stem]; deba@2440: deba@2440: // Finding the node just before the stem node (u) according to deba@2440: // the original thread index deba@2440: for (u = new_stem; thread[u] != stem; u = thread[u]) ; deba@2440: thread[u] = right; deba@2440: deba@2440: // Changing the parent node of stem and shifting stem and deba@2440: // par_stem nodes deba@2440: parent[stem] = par_stem; deba@2440: par_stem = stem; deba@2440: stem = new_stem; deba@2440: deba@2440: // Finding the last successor of stem (u) and the node after it deba@2440: // (right) according to the thread index deba@2440: for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ; deba@2440: right = thread[u]; deba@2440: } deba@2440: parent[u_out] = par_stem; deba@2440: thread[u] = last; deba@2440: deba@2440: if (join == v_out && par_first) { deba@2440: if (first != v_in) thread[first] = right; deba@2440: } else { deba@2440: for (u = v_out; thread[u] != u_out; u = thread[u]) ; deba@2440: thread[u] = right; deba@2440: } deba@2440: } deba@2440: deba@2440: /// \brief Updates pred_edge and forward node maps. deba@2440: void updatePredEdge() { deba@2440: Node u = u_out, v; deba@2440: while (u != u_in) { deba@2440: v = parent[u]; deba@2440: pred_edge[u] = pred_edge[v]; deba@2440: forward[u] = !forward[v]; deba@2440: u = v; deba@2440: } deba@2440: pred_edge[u_in] = in_edge; deba@2440: forward[u_in] = (u_in == graph.source(in_edge)); deba@2440: } deba@2440: deba@2440: /// \brief Updates depth and potential node maps. deba@2440: void updateDepthPotential() { deba@2440: depth[u_in] = depth[v_in] + 1; deba@2440: potential[u_in] = forward[u_in] ? deba@2440: potential[v_in] + cost[pred_edge[u_in]] : deba@2440: potential[v_in] - cost[pred_edge[u_in]]; deba@2440: deba@2440: Node u = thread[u_in], v; deba@2440: while (true) { deba@2440: v = parent[u]; deba@2440: if (v == INVALID) break; deba@2440: depth[u] = depth[v] + 1; deba@2440: potential[u] = forward[u] ? deba@2440: potential[v] + cost[pred_edge[u]] : deba@2440: potential[v] - cost[pred_edge[u]]; deba@2440: if (depth[u] <= depth[v_in]) break; deba@2440: u = thread[u]; deba@2440: } deba@2440: } deba@2440: deba@2440: /// \brief Executes the algorithm. deba@2440: bool start() { deba@2440: // Processing pivots deba@2440: #ifdef _DEBUG_ITER_ deba@2440: int iter_num = 0; deba@2440: #endif deba@2440: #if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT) deba@2440: EdgeIt next_edge(graph); deba@2440: while (findEnteringEdge(next_edge)) deba@2440: #else deba@2440: while (findEnteringEdge()) deba@2440: #endif deba@2440: { deba@2440: join = findJoinNode(); deba@2440: bool change = findLeavingEdge(); deba@2440: changeFlows(change); deba@2440: if (change) { deba@2440: updateThreadParent(); deba@2440: updatePredEdge(); deba@2440: updateDepthPotential(); deba@2440: } deba@2440: #ifdef _DEBUG_ITER_ deba@2440: ++iter_num; deba@2440: #endif deba@2440: } deba@2440: deba@2440: #ifdef _DEBUG_ITER_ deba@2440: std::cout << "Network Simplex algorithm finished. " << iter_num deba@2444: << " pivot iterations performed." << std::endl; deba@2440: #endif deba@2440: deba@2440: // Checking if the flow amount equals zero on all the deba@2440: // artificial edges deba@2440: for (InEdgeIt e(graph, root); e != INVALID; ++e) deba@2440: if (flow[e] > 0) return false; deba@2440: for (OutEdgeIt e(graph, root); e != INVALID; ++e) deba@2440: if (flow[e] > 0) return false; deba@2440: deba@2440: // Copying flow values to flow_result deba@2440: if (lower) { deba@2440: for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) deba@2440: flow_result[e] = (*lower)[e] + flow[edge_ref[e]]; deba@2440: } else { deba@2440: for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) deba@2440: flow_result[e] = flow[edge_ref[e]]; deba@2440: } deba@2440: // Copying potential values to potential_result deba@2440: for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) deba@2440: 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