Revert to long long int since currently I don't know a better solution.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
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.
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
19 #ifndef LEMON_NETWORK_SIMPLEX_H
20 #define LEMON_NETWORK_SIMPLEX_H
22 /// \ingroup min_cost_flow
25 /// \brief Network simplex algorithm for finding a minimum cost flow.
30 #include <lemon/graph_adaptor.h>
31 #include <lemon/graph_utils.h>
32 #include <lemon/smart_graph.h>
33 #include <lemon/math.h>
37 /// \addtogroup min_cost_flow
40 /// \brief Implementation of the network simplex algorithm for
41 /// finding a minimum cost flow.
43 /// \ref NetworkSimplex implements the network simplex algorithm for
44 /// finding a minimum cost flow.
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.
53 /// - Edge capacities and costs should be \e non-negative \e integers.
54 /// - Supply values should be \e signed \e integers.
55 /// - The value types of the maps should be convertible to each other.
56 /// - \c CostMap::Value must be signed type.
58 /// \note \ref NetworkSimplex provides six different pivot rule
59 /// implementations that significantly affect the efficiency of the
61 /// By default a combined pivot rule is used, which is the fastest
62 /// implementation according to our benchmark tests.
63 /// Another pivot rule can be selected using \ref run() function
64 /// with the proper parameter.
66 /// \author Peter Kovacs
68 template < typename Graph,
69 typename LowerMap = typename Graph::template EdgeMap<int>,
70 typename CapacityMap = typename Graph::template EdgeMap<int>,
71 typename CostMap = typename Graph::template EdgeMap<int>,
72 typename SupplyMap = typename Graph::template NodeMap<int> >
75 typedef typename CapacityMap::Value Capacity;
76 typedef typename CostMap::Value Cost;
77 typedef typename SupplyMap::Value Supply;
79 typedef SmartGraph SGraph;
80 GRAPH_TYPEDEFS(typename SGraph);
82 typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
83 typedef typename SGraph::template EdgeMap<Cost> SCostMap;
84 typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
85 typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
87 typedef typename SGraph::template NodeMap<int> IntNodeMap;
88 typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
89 typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
90 typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
91 typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
93 typedef typename Graph::template NodeMap<Node> NodeRefMap;
94 typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
98 /// The type of the flow map.
99 typedef typename Graph::template EdgeMap<Capacity> FlowMap;
100 /// The type of the potential map.
101 typedef typename Graph::template NodeMap<Cost> PotentialMap;
105 /// Enum type to select the pivot rule used by \ref run().
107 FIRST_ELIGIBLE_PIVOT,
110 LIMITED_SEARCH_PIVOT,
111 CANDIDATE_LIST_PIVOT,
117 /// \brief Map adaptor class for handling reduced edge costs.
119 /// Map adaptor class for handling reduced edge costs.
120 class ReducedCostMap : public MapBase<Edge, Cost>
125 const SCostMap &_cost_map;
126 const SPotentialMap &_pot_map;
131 ReducedCostMap( const SGraph &gr,
132 const SCostMap &cost_map,
133 const SPotentialMap &pot_map ) :
134 _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
137 Cost operator[](const Edge &e) const {
138 return _cost_map[e] + _pot_map[_gr.source(e)]
139 - _pot_map[_gr.target(e)];
142 }; //class ReducedCostMap
146 /// \brief Implementation of the "First Eligible" pivot rule for the
147 /// \ref NetworkSimplex "network simplex" algorithm.
149 /// This class implements the "First Eligible" pivot rule
150 /// for the \ref NetworkSimplex "network simplex" algorithm.
151 class FirstEligiblePivotRule
161 FirstEligiblePivotRule(NetworkSimplex &ns) :
162 _ns(ns), _next_edge(ns._graph) {}
164 /// Finds the next entering edge.
165 bool findEnteringEdge() {
166 for (EdgeIt e = _next_edge; e != INVALID; ++e) {
167 if (_ns._state[e] * _ns._red_cost[e] < 0) {
173 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
174 if (_ns._state[e] * _ns._red_cost[e] < 0) {
182 }; //class FirstEligiblePivotRule
184 /// \brief Implementation of the "Best Eligible" pivot rule for the
185 /// \ref NetworkSimplex "network simplex" algorithm.
187 /// This class implements the "Best Eligible" pivot rule
188 /// for the \ref NetworkSimplex "network simplex" algorithm.
189 class BestEligiblePivotRule
198 BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
200 /// Finds the next entering edge.
201 bool findEnteringEdge() {
203 for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
204 if (_ns._state[e] * _ns._red_cost[e] < min) {
205 min = _ns._state[e] * _ns._red_cost[e];
211 }; //class BestEligiblePivotRule
213 /// \brief Implementation of the "Block Search" pivot rule for the
214 /// \ref NetworkSimplex "network simplex" algorithm.
216 /// This class implements the "Block Search" pivot rule
217 /// for the \ref NetworkSimplex "network simplex" algorithm.
218 class BlockSearchPivotRule
223 EdgeIt _next_edge, _min_edge;
226 static const int MIN_BLOCK_SIZE = 10;
231 BlockSearchPivotRule(NetworkSimplex &ns) :
232 _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
234 _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
235 if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
238 /// Finds the next entering edge.
239 bool findEnteringEdge() {
242 for (EdgeIt e = _next_edge; e != INVALID; ++e) {
243 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
247 if (++cnt == _block_size) {
253 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
254 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
258 if (++cnt == _block_size) {
264 _ns._in_edge = _min_edge;
265 _next_edge = ++_min_edge;
268 }; //class BlockSearchPivotRule
270 /// \brief Implementation of the "Limited Search" pivot rule for the
271 /// \ref NetworkSimplex "network simplex" algorithm.
273 /// This class implements the "Limited Search" pivot rule
274 /// for the \ref NetworkSimplex "network simplex" algorithm.
275 class LimitedSearchPivotRule
280 EdgeIt _next_edge, _min_edge;
283 static const int MIN_SAMPLE_SIZE = 10;
284 static const double SAMPLE_SIZE_FACTOR = 0.0015;
289 LimitedSearchPivotRule(NetworkSimplex &ns) :
290 _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
292 _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph));
293 if (_sample_size < MIN_SAMPLE_SIZE)
294 _sample_size = MIN_SAMPLE_SIZE;
297 /// Finds the next entering edge.
298 bool findEnteringEdge() {
301 for (EdgeIt e = _next_edge; e != INVALID; ++e) {
302 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
306 if (curr < 0 && ++cnt == _sample_size) break;
309 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
310 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
314 if (curr < 0 && ++cnt == _sample_size) break;
317 _ns._in_edge = _min_edge;
318 _next_edge = ++_min_edge;
321 }; //class LimitedSearchPivotRule
323 /// \brief Implementation of the "Candidate List" pivot rule for the
324 /// \ref NetworkSimplex "network simplex" algorithm.
326 /// This class implements the "Candidate List" pivot rule
327 /// for the \ref NetworkSimplex "network simplex" algorithm.
328 class CandidateListPivotRule
334 // The list of candidate edges.
335 std::vector<Edge> _candidates;
336 // The maximum length of the edge list.
338 // The maximum number of minor iterations between two major
345 static const double LIST_LENGTH_FACTOR = 0.002;
346 static const double MINOR_LIMIT_FACTOR = 0.1;
347 static const int MIN_LIST_LENGTH = 10;
348 static const int MIN_MINOR_LIMIT = 2;
353 CandidateListPivotRule(NetworkSimplex &ns) :
354 _ns(ns), _next_edge(ns._graph)
356 int edge_num = countEdges(_ns._graph);
358 _list_length = int(edge_num * LIST_LENGTH_FACTOR);
359 if (_list_length < MIN_LIST_LENGTH)
360 _list_length = MIN_LIST_LENGTH;
361 _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR);
362 if (_minor_limit < MIN_MINOR_LIMIT)
363 _minor_limit = MIN_MINOR_LIMIT;
366 /// Finds the next entering edge.
367 bool findEnteringEdge() {
369 if (_minor_count < _minor_limit && _candidates.size() > 0) {
374 for (int i = 0; i < int(_candidates.size()); ++i) {
376 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
381 if (min < 0) return true;
386 EdgeIt e = _next_edge;
388 for ( ; e != INVALID; ++e) {
389 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
390 _candidates.push_back(e);
395 if (int(_candidates.size()) == _list_length) break;
398 if (int(_candidates.size()) < _list_length) {
399 for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
400 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
401 _candidates.push_back(e);
406 if (int(_candidates.size()) == _list_length) break;
410 if (_candidates.size() == 0) return false;
415 }; //class CandidateListPivotRule
419 // State constants for edges
426 // Constant for the combined pivot rule.
427 static const int COMBINED_PIVOT_MAX_DEG = 5;
431 // The directed graph the algorithm runs on
433 // The original graph
434 const Graph &_graph_ref;
435 // The original lower bound map
436 const LowerMap *_lower;
438 SCapacityMap _capacity;
445 // Edge map of the current flow
447 // Node map of the current potentials
448 SPotentialMap _potential;
450 // The depth node map of the spanning tree structure
452 // The parent node map of the spanning tree structure
454 // The pred_edge node map of the spanning tree structure
455 EdgeNodeMap _pred_edge;
456 // The thread node map of the spanning tree structure
458 // The forward node map of the spanning tree structure
459 BoolNodeMap _forward;
460 // The state edge map
462 // The root node of the starting spanning tree
465 // The reduced cost map
466 ReducedCostMap _red_cost;
468 // Members for handling the original graph
469 FlowMap *_flow_result;
470 PotentialMap *_potential_result;
472 bool _local_potential;
473 NodeRefMap _node_ref;
474 EdgeRefMap _edge_ref;
476 // The entering edge of the current pivot iteration.
479 // Temporary nodes used in the current pivot iteration.
480 Node join, u_in, v_in, u_out, v_out;
481 Node right, first, second, last;
482 Node stem, par_stem, new_stem;
483 // The maximum augment amount along the found cycle in the current
489 /// \brief General constructor (with lower bounds).
491 /// General constructor (with lower bounds).
493 /// \param graph The directed graph the algorithm runs on.
494 /// \param lower The lower bounds of the edges.
495 /// \param capacity The capacities (upper bounds) of the edges.
496 /// \param cost The cost (length) values of the edges.
497 /// \param supply The supply values of the nodes (signed).
498 NetworkSimplex( const Graph &graph,
499 const LowerMap &lower,
500 const CapacityMap &capacity,
502 const SupplyMap &supply ) :
503 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
504 _cost(_graph), _supply(_graph), _flow(_graph),
505 _potential(_graph), _depth(_graph), _parent(_graph),
506 _pred_edge(_graph), _thread(_graph), _forward(_graph),
507 _state(_graph), _red_cost(_graph, _cost, _potential),
508 _flow_result(0), _potential_result(0),
509 _local_flow(false), _local_potential(false),
510 _node_ref(graph), _edge_ref(graph)
512 // Checking the sum of supply values
514 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
516 if (!(_valid_supply = sum == 0)) return;
518 // Copying _graph_ref to _graph
519 _graph.reserveNode(countNodes(_graph_ref) + 1);
520 _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
521 copyGraph(_graph, _graph_ref)
522 .edgeMap(_cost, cost)
527 // Removing non-zero lower bounds
528 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
529 _capacity[_edge_ref[e]] = capacity[e] - lower[e];
531 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
532 Supply s = supply[n];
533 for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
535 for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
537 _supply[_node_ref[n]] = s;
541 /// \brief General constructor (without lower bounds).
543 /// General constructor (without lower bounds).
545 /// \param graph The directed graph the algorithm runs on.
546 /// \param capacity The capacities (upper bounds) of the edges.
547 /// \param cost The cost (length) values of the edges.
548 /// \param supply The supply values of the nodes (signed).
549 NetworkSimplex( const Graph &graph,
550 const CapacityMap &capacity,
552 const SupplyMap &supply ) :
553 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
554 _cost(_graph), _supply(_graph), _flow(_graph),
555 _potential(_graph), _depth(_graph), _parent(_graph),
556 _pred_edge(_graph), _thread(_graph), _forward(_graph),
557 _state(_graph), _red_cost(_graph, _cost, _potential),
558 _flow_result(0), _potential_result(0),
559 _local_flow(false), _local_potential(false),
560 _node_ref(graph), _edge_ref(graph)
562 // Checking the sum of supply values
564 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
566 if (!(_valid_supply = sum == 0)) return;
568 // Copying _graph_ref to graph
569 copyGraph(_graph, _graph_ref)
570 .edgeMap(_capacity, capacity)
571 .edgeMap(_cost, cost)
572 .nodeMap(_supply, supply)
578 /// \brief Simple constructor (with lower bounds).
580 /// Simple constructor (with lower bounds).
582 /// \param graph The directed graph the algorithm runs on.
583 /// \param lower The lower bounds of the edges.
584 /// \param capacity The capacities (upper bounds) of the edges.
585 /// \param cost The cost (length) values of the edges.
586 /// \param s The source node.
587 /// \param t The target node.
588 /// \param flow_value The required amount of flow from node \c s
589 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
590 NetworkSimplex( const Graph &graph,
591 const LowerMap &lower,
592 const CapacityMap &capacity,
594 typename Graph::Node s,
595 typename Graph::Node t,
596 typename SupplyMap::Value flow_value ) :
597 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
598 _cost(_graph), _supply(_graph), _flow(_graph),
599 _potential(_graph), _depth(_graph), _parent(_graph),
600 _pred_edge(_graph), _thread(_graph), _forward(_graph),
601 _state(_graph), _red_cost(_graph, _cost, _potential),
602 _flow_result(0), _potential_result(0),
603 _local_flow(false), _local_potential(false),
604 _node_ref(graph), _edge_ref(graph)
606 // Copying _graph_ref to graph
607 copyGraph(_graph, _graph_ref)
608 .edgeMap(_cost, cost)
613 // Removing non-zero lower bounds
614 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
615 _capacity[_edge_ref[e]] = capacity[e] - lower[e];
617 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
619 if (n == s) sum = flow_value;
620 if (n == t) sum = -flow_value;
621 for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
623 for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
625 _supply[_node_ref[n]] = sum;
627 _valid_supply = true;
630 /// \brief Simple constructor (without lower bounds).
632 /// Simple constructor (without lower bounds).
634 /// \param graph The directed graph the algorithm runs on.
635 /// \param capacity The capacities (upper bounds) of the edges.
636 /// \param cost The cost (length) values of the edges.
637 /// \param s The source node.
638 /// \param t The target node.
639 /// \param flow_value The required amount of flow from node \c s
640 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
641 NetworkSimplex( const Graph &graph,
642 const CapacityMap &capacity,
644 typename Graph::Node s,
645 typename Graph::Node t,
646 typename SupplyMap::Value flow_value ) :
647 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
648 _cost(_graph), _supply(_graph, 0), _flow(_graph),
649 _potential(_graph), _depth(_graph), _parent(_graph),
650 _pred_edge(_graph), _thread(_graph), _forward(_graph),
651 _state(_graph), _red_cost(_graph, _cost, _potential),
652 _flow_result(0), _potential_result(0),
653 _local_flow(false), _local_potential(false),
654 _node_ref(graph), _edge_ref(graph)
656 // Copying _graph_ref to graph
657 copyGraph(_graph, _graph_ref)
658 .edgeMap(_capacity, capacity)
659 .edgeMap(_cost, cost)
663 _supply[_node_ref[s]] = flow_value;
664 _supply[_node_ref[t]] = -flow_value;
665 _valid_supply = true;
670 if (_local_flow) delete _flow_result;
671 if (_local_potential) delete _potential_result;
674 /// \brief Sets the flow map.
676 /// Sets the flow map.
678 /// \return \c (*this)
679 NetworkSimplex& flowMap(FlowMap &map) {
688 /// \brief Sets the potential map.
690 /// Sets the potential map.
692 /// \return \c (*this)
693 NetworkSimplex& potentialMap(PotentialMap &map) {
694 if (_local_potential) {
695 delete _potential_result;
696 _local_potential = false;
698 _potential_result = ↦
702 /// \name Execution control
703 /// The only way to execute the algorithm is to call the run()
708 /// \brief Runs the algorithm.
710 /// Runs the algorithm.
712 /// \param pivot_rule The pivot rule that is used during the
715 /// The available pivot rules:
717 /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
718 /// a wraparound fashion in every iteration
719 /// (\ref FirstEligiblePivotRule).
721 /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
722 /// every iteration (\ref BestEligiblePivotRule).
724 /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
725 /// every iteration in a wraparound fashion and the best eligible
726 /// edge is selected from this block (\ref BlockSearchPivotRule).
728 /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
729 /// examined in every iteration in a wraparound fashion and the best
730 /// one is selected from them (\ref LimitedSearchPivotRule).
732 /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
733 /// built from eligible edges and it is used for edge selection in
734 /// the following minor iterations (\ref CandidateListPivotRule).
736 /// - COMBINED_PIVOT This is a combined version of the two fastest
738 /// For rather sparse graphs \ref LimitedSearchPivotRule
739 /// "Limited Search" implementation is used, otherwise
740 /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
741 /// According to our benchmark tests this combined method is the
744 /// \return \c true if a feasible flow can be found.
745 bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
746 return init() && start(pivot_rule);
751 /// \name Query Functions
752 /// The result of the algorithm can be obtained using these
754 /// \n run() must be called before using them.
758 /// \brief Returns a const reference to the edge map storing the
761 /// Returns a const reference to the edge map storing the found flow.
763 /// \pre \ref run() must be called before using this function.
764 const FlowMap& flowMap() const {
765 return *_flow_result;
768 /// \brief Returns a const reference to the node map storing the
769 /// found potentials (the dual solution).
771 /// Returns a const reference to the node map storing the found
772 /// potentials (the dual solution).
774 /// \pre \ref run() must be called before using this function.
775 const PotentialMap& potentialMap() const {
776 return *_potential_result;
779 /// \brief Returns the flow on the given edge.
781 /// Returns the flow on the given edge.
783 /// \pre \ref run() must be called before using this function.
784 Capacity flow(const typename Graph::Edge& edge) const {
785 return (*_flow_result)[edge];
788 /// \brief Returns the potential of the given node.
790 /// Returns the potential of the given node.
792 /// \pre \ref run() must be called before using this function.
793 Cost potential(const typename Graph::Node& node) const {
794 return (*_potential_result)[node];
797 /// \brief Returns the total cost of the found flow.
799 /// Returns the total cost of the found flow. The complexity of the
800 /// function is \f$ O(e) \f$.
802 /// \pre \ref run() must be called before using this function.
803 Cost totalCost() const {
805 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
806 c += (*_flow_result)[e] * _cost[_edge_ref[e]];
814 /// \brief Extends the underlaying graph and initializes all the
815 /// node and edge maps.
817 if (!_valid_supply) return false;
819 // Initializing result maps
821 _flow_result = new FlowMap(_graph_ref);
824 if (!_potential_result) {
825 _potential_result = new PotentialMap(_graph_ref);
826 _local_potential = true;
829 // Initializing state and flow maps
830 for (EdgeIt e(_graph); e != INVALID; ++e) {
832 _state[e] = STATE_LOWER;
835 // Adding an artificial root node to the graph
836 _root = _graph.addNode();
837 _parent[_root] = INVALID;
838 _pred_edge[_root] = INVALID;
841 _potential[_root] = 0;
843 // Adding artificial edges to the graph and initializing the node
844 // maps of the spanning tree data structure
847 Cost max_cost = std::numeric_limits<Cost>::max() / 4;
848 for (NodeIt u(_graph); u != INVALID; ++u) {
849 if (u == _root) continue;
854 if (_supply[u] >= 0) {
855 e = _graph.addEdge(u, _root);
856 _flow[e] = _supply[u];
858 _potential[u] = -max_cost;
860 e = _graph.addEdge(_root, u);
861 _flow[e] = -_supply[u];
863 _potential[u] = max_cost;
866 _capacity[e] = std::numeric_limits<Capacity>::max();
867 _state[e] = STATE_TREE;
870 _thread[last] = _root;
875 /// Finds the join node.
876 Node findJoinNode() {
877 Node u = _graph.source(_in_edge);
878 Node v = _graph.target(_in_edge);
880 if (_depth[u] == _depth[v]) {
884 else if (_depth[u] > _depth[v]) u = _parent[u];
890 /// \brief Finds the leaving edge of the cycle. Returns \c true if
891 /// the leaving edge is not the same as the entering edge.
892 bool findLeavingEdge() {
893 // Initializing first and second nodes according to the direction
895 if (_state[_in_edge] == STATE_LOWER) {
896 first = _graph.source(_in_edge);
897 second = _graph.target(_in_edge);
899 first = _graph.target(_in_edge);
900 second = _graph.source(_in_edge);
902 delta = _capacity[_in_edge];
907 // Searching the cycle along the path form the first node to the
909 for (Node u = first; u != join; u = _parent[u]) {
911 d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
920 // Searching the cycle along the path form the second node to the
922 for (Node u = second; u != join; u = _parent[u]) {
924 d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
936 /// Changes \c flow and \c state edge maps.
937 void changeFlows(bool change) {
938 // Augmenting along the cycle
940 Capacity val = _state[_in_edge] * delta;
941 _flow[_in_edge] += val;
942 for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
943 _flow[_pred_edge[u]] += _forward[u] ? -val : val;
945 for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
946 _flow[_pred_edge[u]] += _forward[u] ? val : -val;
949 // Updating the state of the entering and leaving edges
951 _state[_in_edge] = STATE_TREE;
952 _state[_pred_edge[u_out]] =
953 (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
955 _state[_in_edge] = -_state[_in_edge];
959 /// Updates \c thread and \c parent node maps.
960 void updateThreadParent() {
962 v_out = _parent[u_out];
964 // Handling the case when join and v_out coincide
965 bool par_first = false;
967 for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
970 while (_thread[u] != u_out) u = _thread[u];
975 // Finding the last successor of u_in (u) and the node after it
976 // (right) according to the thread index
977 for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
979 if (_thread[v_in] == u_out) {
980 for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
981 if (last == u_out) last = _thread[last];
983 else last = _thread[v_in];
985 // Updating stem nodes
986 _thread[v_in] = stem = u_in;
988 while (stem != u_out) {
989 _thread[u] = new_stem = _parent[stem];
991 // Finding the node just before the stem node (u) according to
992 // the original thread index
993 for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
996 // Changing the parent node of stem and shifting stem and
998 _parent[stem] = par_stem;
1002 // Finding the last successor of stem (u) and the node after it
1003 // (right) according to the thread index
1004 for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
1007 _parent[u_out] = par_stem;
1010 if (join == v_out && par_first) {
1011 if (first != v_in) _thread[first] = right;
1013 for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
1018 /// Updates \c pred_edge and \c forward node maps.
1019 void updatePredEdge() {
1023 _pred_edge[u] = _pred_edge[v];
1024 _forward[u] = !_forward[v];
1027 _pred_edge[u_in] = _in_edge;
1028 _forward[u_in] = (u_in == _graph.source(_in_edge));
1031 /// Updates \c depth and \c potential node maps.
1032 void updateDepthPotential() {
1033 _depth[u_in] = _depth[v_in] + 1;
1034 _potential[u_in] = _forward[u_in] ?
1035 _potential[v_in] - _cost[_pred_edge[u_in]] :
1036 _potential[v_in] + _cost[_pred_edge[u_in]];
1038 Node u = _thread[u_in], v;
1041 if (v == INVALID) break;
1042 _depth[u] = _depth[v] + 1;
1043 _potential[u] = _forward[u] ?
1044 _potential[v] - _cost[_pred_edge[u]] :
1045 _potential[v] + _cost[_pred_edge[u]];
1046 if (_depth[u] <= _depth[v_in]) break;
1051 /// Executes the algorithm.
1052 bool start(PivotRuleEnum pivot_rule) {
1053 switch (pivot_rule) {
1054 case FIRST_ELIGIBLE_PIVOT:
1055 return start<FirstEligiblePivotRule>();
1056 case BEST_ELIGIBLE_PIVOT:
1057 return start<BestEligiblePivotRule>();
1058 case BLOCK_SEARCH_PIVOT:
1059 return start<BlockSearchPivotRule>();
1060 case LIMITED_SEARCH_PIVOT:
1061 return start<LimitedSearchPivotRule>();
1062 case CANDIDATE_LIST_PIVOT:
1063 return start<CandidateListPivotRule>();
1064 case COMBINED_PIVOT:
1065 if ( countEdges(_graph) / countNodes(_graph) <=
1066 COMBINED_PIVOT_MAX_DEG )
1067 return start<LimitedSearchPivotRule>();
1069 return start<BlockSearchPivotRule>();
1074 template<class PivotRuleImplementation>
1076 PivotRuleImplementation pivot(*this);
1078 // Executing the network simplex algorithm
1079 while (pivot.findEnteringEdge()) {
1080 join = findJoinNode();
1081 bool change = findLeavingEdge();
1082 changeFlows(change);
1084 updateThreadParent();
1086 updateDepthPotential();
1090 // Checking if the flow amount equals zero on all the artificial
1092 for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
1093 if (_flow[e] > 0) return false;
1094 for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
1095 if (_flow[e] > 0) return false;
1097 // Copying flow values to _flow_result
1099 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1100 (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
1102 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1103 (*_flow_result)[e] = _flow[_edge_ref[e]];
1105 // Copying potential values to _potential_result
1106 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1107 (*_potential_result)[n] = _potential[_node_ref[n]];
1112 }; //class NetworkSimplex
1118 #endif //LEMON_NETWORK_SIMPLEX_H