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.
31 #include <lemon/graph_utils.h>
32 #include <lemon/math.h>
36 /// \addtogroup min_cost_flow
39 /// \brief Implementation of the primal network simplex algorithm
40 /// for finding a minimum cost flow.
42 /// \ref NetworkSimplex implements the primal network simplex algorithm
43 /// for finding a minimum cost flow.
45 /// \tparam Graph The directed graph type the algorithm runs on.
46 /// \tparam LowerMap The type of the lower bound map.
47 /// \tparam CapacityMap The type of the capacity (upper bound) map.
48 /// \tparam CostMap The type of the cost (length) map.
49 /// \tparam SupplyMap The type of the supply map.
52 /// - Edge capacities and costs should be \e non-negative \e integers.
53 /// - Supply values should be \e signed \e integers.
54 /// - The value types of the maps should be convertible to each other.
55 /// - \c CostMap::Value must be signed type.
57 /// \note \ref NetworkSimplex provides five different pivot rule
58 /// implementations that significantly affect the efficiency of the
60 /// By default "Block Search" pivot rule is used, which proved to be
61 /// by far the most efficient according to our benchmark tests.
62 /// However another pivot rule can be selected using \ref run()
63 /// function with the proper parameter.
65 /// \author Peter Kovacs
66 template < typename Graph,
67 typename LowerMap = typename Graph::template EdgeMap<int>,
68 typename CapacityMap = typename Graph::template EdgeMap<int>,
69 typename CostMap = typename Graph::template EdgeMap<int>,
70 typename SupplyMap = typename Graph::template NodeMap<int> >
73 GRAPH_TYPEDEFS(typename Graph);
75 typedef typename CapacityMap::Value Capacity;
76 typedef typename CostMap::Value Cost;
77 typedef typename SupplyMap::Value Supply;
79 typedef std::vector<Edge> EdgeVector;
80 typedef std::vector<Node> NodeVector;
81 typedef std::vector<int> IntVector;
82 typedef std::vector<bool> BoolVector;
83 typedef std::vector<Capacity> CapacityVector;
84 typedef std::vector<Cost> CostVector;
85 typedef std::vector<Supply> SupplyVector;
87 typedef typename Graph::template NodeMap<int> IntNodeMap;
91 /// The type of the flow map.
92 typedef typename Graph::template EdgeMap<Capacity> FlowMap;
93 /// The type of the potential map.
94 typedef typename Graph::template NodeMap<Cost> PotentialMap;
98 /// Enum type to select the pivot rule used by \ref run().
100 FIRST_ELIGIBLE_PIVOT,
103 CANDIDATE_LIST_PIVOT,
109 /// \brief Implementation of the "First Eligible" pivot rule for the
110 /// \ref NetworkSimplex "network simplex" algorithm.
112 /// This class implements the "First Eligible" pivot rule
113 /// for the \ref NetworkSimplex "network simplex" algorithm.
115 /// For more information see \ref NetworkSimplex::run().
116 class FirstEligiblePivotRule
120 // References to the NetworkSimplex class
121 const EdgeVector &_edge;
122 const IntVector &_source;
123 const IntVector &_target;
124 const CostVector &_cost;
125 const IntVector &_state;
126 const CostVector &_pi;
136 FirstEligiblePivotRule(NetworkSimplex &ns) :
137 _edge(ns._edge), _source(ns._source), _target(ns._target),
138 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
139 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
142 /// Find next entering edge
143 bool findEnteringEdge() {
145 for (int e = _next_edge; e < _edge_num; ++e) {
146 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
153 for (int e = 0; e < _next_edge; ++e) {
154 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
164 }; //class FirstEligiblePivotRule
167 /// \brief Implementation of the "Best Eligible" pivot rule for the
168 /// \ref NetworkSimplex "network simplex" algorithm.
170 /// This class implements the "Best Eligible" pivot rule
171 /// for the \ref NetworkSimplex "network simplex" algorithm.
173 /// For more information see \ref NetworkSimplex::run().
174 class BestEligiblePivotRule
178 // References to the NetworkSimplex class
179 const EdgeVector &_edge;
180 const IntVector &_source;
181 const IntVector &_target;
182 const CostVector &_cost;
183 const IntVector &_state;
184 const CostVector &_pi;
191 BestEligiblePivotRule(NetworkSimplex &ns) :
192 _edge(ns._edge), _source(ns._source), _target(ns._target),
193 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
194 _in_edge(ns._in_edge), _edge_num(ns._edge_num)
197 /// Find next entering edge
198 bool findEnteringEdge() {
200 for (int e = 0; e < _edge_num; ++e) {
201 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
210 }; //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.
219 /// For more information see \ref NetworkSimplex::run().
220 class BlockSearchPivotRule
224 // References to the NetworkSimplex class
225 const EdgeVector &_edge;
226 const IntVector &_source;
227 const IntVector &_target;
228 const CostVector &_cost;
229 const IntVector &_state;
230 const CostVector &_pi;
241 BlockSearchPivotRule(NetworkSimplex &ns) :
242 _edge(ns._edge), _source(ns._source), _target(ns._target),
243 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
244 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
246 // The main parameters of the pivot rule
247 const double BLOCK_SIZE_FACTOR = 0.5;
248 const int MIN_BLOCK_SIZE = 10;
250 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
254 /// Find next entering edge
255 bool findEnteringEdge() {
257 int cnt = _block_size;
259 for (e = _next_edge; e < _edge_num; ++e) {
260 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
266 if (min < 0) goto search_end;
270 for (e = 0; e < _next_edge; ++e) {
271 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
277 if (min < 0) goto search_end;
281 if (min >= 0) return false;
288 }; //class BlockSearchPivotRule
291 /// \brief Implementation of the "Candidate List" pivot rule for the
292 /// \ref NetworkSimplex "network simplex" algorithm.
294 /// This class implements the "Candidate List" pivot rule
295 /// for the \ref NetworkSimplex "network simplex" algorithm.
297 /// For more information see \ref NetworkSimplex::run().
298 class CandidateListPivotRule
302 // References to the NetworkSimplex class
303 const EdgeVector &_edge;
304 const IntVector &_source;
305 const IntVector &_target;
306 const CostVector &_cost;
307 const IntVector &_state;
308 const CostVector &_pi;
313 IntVector _candidates;
314 int _list_length, _minor_limit;
315 int _curr_length, _minor_count;
321 CandidateListPivotRule(NetworkSimplex &ns) :
322 _edge(ns._edge), _source(ns._source), _target(ns._target),
323 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
324 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
326 // The main parameters of the pivot rule
327 const double LIST_LENGTH_FACTOR = 0.25;
328 const int MIN_LIST_LENGTH = 10;
329 const double MINOR_LIMIT_FACTOR = 0.1;
330 const int MIN_MINOR_LIMIT = 3;
332 _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)),
334 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
336 _curr_length = _minor_count = 0;
337 _candidates.resize(_list_length);
340 /// Find next entering edge
341 bool findEnteringEdge() {
344 if (_curr_length > 0 && _minor_count < _minor_limit) {
345 // Minor iteration: select the best eligible edge from the
346 // current candidate list
349 for (int i = 0; i < _curr_length; ++i) {
351 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
357 _candidates[i--] = _candidates[--_curr_length];
360 if (min < 0) return true;
363 // Major iteration: build a new candidate list
366 for (e = _next_edge; e < _edge_num; ++e) {
367 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
369 _candidates[_curr_length++] = e;
374 if (_curr_length == _list_length) goto search_end;
377 for (e = 0; e < _next_edge; ++e) {
378 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
380 _candidates[_curr_length++] = e;
385 if (_curr_length == _list_length) goto search_end;
388 if (_curr_length == 0) return false;
396 }; //class CandidateListPivotRule
399 /// \brief Implementation of the "Altering Candidate List" pivot rule
400 /// for the \ref NetworkSimplex "network simplex" algorithm.
402 /// This class implements the "Altering Candidate List" pivot rule
403 /// for the \ref NetworkSimplex "network simplex" algorithm.
405 /// For more information see \ref NetworkSimplex::run().
406 class AlteringListPivotRule
410 // References to the NetworkSimplex class
411 const EdgeVector &_edge;
412 const IntVector &_source;
413 const IntVector &_target;
414 const CostVector &_cost;
415 const IntVector &_state;
416 const CostVector &_pi;
420 int _block_size, _head_length, _curr_length;
422 IntVector _candidates;
423 CostVector _cand_cost;
425 // Functor class to compare edges during sort of the candidate list
429 const CostVector &_map;
431 SortFunc(const CostVector &map) : _map(map) {}
432 bool operator()(int left, int right) {
433 return _map[left] > _map[right];
442 AlteringListPivotRule(NetworkSimplex &ns) :
443 _edge(ns._edge), _source(ns._source), _target(ns._target),
444 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
445 _in_edge(ns._in_edge), _edge_num(ns._edge_num),
446 _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
448 // The main parameters of the pivot rule
449 const double BLOCK_SIZE_FACTOR = 1.0;
450 const int MIN_BLOCK_SIZE = 10;
451 const double HEAD_LENGTH_FACTOR = 0.1;
452 const int MIN_HEAD_LENGTH = 3;
454 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
456 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
458 _candidates.resize(_head_length + _block_size);
462 /// Find next entering edge
463 bool findEnteringEdge() {
464 // Check the current candidate list
466 for (int i = 0; i < _curr_length; ++i) {
468 _cand_cost[e] = _state[e] *
469 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
470 if (_cand_cost[e] >= 0) {
471 _candidates[i--] = _candidates[--_curr_length];
476 int cnt = _block_size;
477 int limit = _head_length;
479 for (e = _next_edge; e < _edge_num; ++e) {
480 _cand_cost[e] = _state[e] *
481 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
482 if (_cand_cost[e] < 0) {
483 _candidates[_curr_length++] = e;
486 if (_curr_length > limit) goto search_end;
491 for (e = 0; e < _next_edge; ++e) {
492 _cand_cost[e] = _state[e] *
493 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
494 if (_cand_cost[e] < 0) {
495 _candidates[_curr_length++] = e;
498 if (_curr_length > limit) goto search_end;
503 if (_curr_length == 0) return false;
507 // Make heap of the candidate list (approximating a partial sort)
508 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
511 // Pop the first element of the heap
512 _in_edge = _candidates[0];
514 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
516 _curr_length = std::min(_head_length, _curr_length - 1);
520 }; //class AlteringListPivotRule
524 // State constants for edges
533 // The original graph
534 const Graph &_orig_graph;
535 // The original lower bound map
536 const LowerMap *_orig_lower;
537 // The original capacity map
538 const CapacityMap &_orig_cap;
539 // The original cost map
540 const CostMap &_orig_cost;
541 // The original supply map
542 const SupplyMap *_orig_supply;
543 // The source node (if no supply map was given)
545 // The target node (if no supply map was given)
547 // The flow value (if no supply map was given)
548 Capacity _orig_flow_value;
550 // The flow result map
551 FlowMap *_flow_result;
552 // The potential result map
553 PotentialMap *_potential_result;
554 // Indicate if the flow result map is local
556 // Indicate if the potential result map is local
557 bool _local_potential;
559 // The edge references
561 // The node references
565 // The source nodes of the edges
567 // The target nodess of the edges
570 // The (modified) capacity vector
574 // The (modified) supply vector
576 // The current flow vector
577 CapacityVector _flow;
578 // The current potential vector
581 // The number of nodes in the original graph
583 // The number of edges in the original graph
586 // The parent vector of the spanning tree structure
588 // The pred_edge vector of the spanning tree structure
590 // The thread vector of the spanning tree structure
593 IntVector _rev_thread;
595 IntVector _last_succ;
597 IntVector _dirty_revs;
599 // The forward vector of the spanning tree structure
606 // The entering edge of the current pivot iteration
609 // Temporary nodes used in the current pivot iteration
610 int join, u_in, v_in, u_out, v_out;
611 int right, first, second, last;
612 int stem, par_stem, new_stem;
614 // The maximum augment amount along the found cycle in the current
620 /// \brief General constructor (with lower bounds).
622 /// General constructor (with lower bounds).
624 /// \param graph The directed graph the algorithm runs on.
625 /// \param lower The lower bounds of the edges.
626 /// \param capacity The capacities (upper bounds) of the edges.
627 /// \param cost The cost (length) values of the edges.
628 /// \param supply The supply values of the nodes (signed).
629 NetworkSimplex( const Graph &graph,
630 const LowerMap &lower,
631 const CapacityMap &capacity,
633 const SupplyMap &supply ) :
634 _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
635 _orig_cost(cost), _orig_supply(&supply),
636 _flow_result(NULL), _potential_result(NULL),
637 _local_flow(false), _local_potential(false),
641 /// \brief General constructor (without lower bounds).
643 /// General constructor (without lower bounds).
645 /// \param graph The directed graph the algorithm runs on.
646 /// \param capacity The capacities (upper bounds) of the edges.
647 /// \param cost The cost (length) values of the edges.
648 /// \param supply The supply values of the nodes (signed).
649 NetworkSimplex( const Graph &graph,
650 const CapacityMap &capacity,
652 const SupplyMap &supply ) :
653 _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
654 _orig_cost(cost), _orig_supply(&supply),
655 _flow_result(NULL), _potential_result(NULL),
656 _local_flow(false), _local_potential(false),
660 /// \brief Simple constructor (with lower bounds).
662 /// Simple constructor (with lower bounds).
664 /// \param graph The directed graph the algorithm runs on.
665 /// \param lower The lower bounds of the edges.
666 /// \param capacity The capacities (upper bounds) of the edges.
667 /// \param cost The cost (length) values of the edges.
668 /// \param s The source node.
669 /// \param t The target node.
670 /// \param flow_value The required amount of flow from node \c s
671 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
672 NetworkSimplex( const Graph &graph,
673 const LowerMap &lower,
674 const CapacityMap &capacity,
677 Capacity flow_value ) :
678 _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
679 _orig_cost(cost), _orig_supply(NULL),
680 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
681 _flow_result(NULL), _potential_result(NULL),
682 _local_flow(false), _local_potential(false),
686 /// \brief Simple constructor (without lower bounds).
688 /// Simple constructor (without lower bounds).
690 /// \param graph The directed graph the algorithm runs on.
691 /// \param capacity The capacities (upper bounds) of the edges.
692 /// \param cost The cost (length) values of the edges.
693 /// \param s The source node.
694 /// \param t The target node.
695 /// \param flow_value The required amount of flow from node \c s
696 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
697 NetworkSimplex( const Graph &graph,
698 const CapacityMap &capacity,
701 Capacity flow_value ) :
702 _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
703 _orig_cost(cost), _orig_supply(NULL),
704 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
705 _flow_result(NULL), _potential_result(NULL),
706 _local_flow(false), _local_potential(false),
712 if (_local_flow) delete _flow_result;
713 if (_local_potential) delete _potential_result;
716 /// \brief Set the flow map.
718 /// Set the flow map.
720 /// \return \c (*this)
721 NetworkSimplex& flowMap(FlowMap &map) {
730 /// \brief Set the potential map.
732 /// Set the potential map.
734 /// \return \c (*this)
735 NetworkSimplex& potentialMap(PotentialMap &map) {
736 if (_local_potential) {
737 delete _potential_result;
738 _local_potential = false;
740 _potential_result = ↦
744 /// \name Execution control
748 /// \brief Runs the algorithm.
750 /// Runs the algorithm.
752 /// \param pivot_rule The pivot rule that is used during the
755 /// The available pivot rules:
757 /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
758 /// a wraparound fashion in every iteration
759 /// (\ref FirstEligiblePivotRule).
761 /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
762 /// every iteration (\ref BestEligiblePivotRule).
764 /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
765 /// every iteration in a wraparound fashion and the best eligible
766 /// edge is selected from this block (\ref BlockSearchPivotRule).
768 /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
769 /// built from eligible edges in a wraparound fashion and in the
770 /// following minor iterations the best eligible edge is selected
771 /// from this list (\ref CandidateListPivotRule).
773 /// - ALTERING_LIST_PIVOT It is a modified version of the
774 /// "Candidate List" pivot rule. It keeps only the several best
775 /// eligible edges from the former candidate list and extends this
776 /// list in every iteration (\ref AlteringListPivotRule).
778 /// According to our comprehensive benchmark tests the "Block Search"
779 /// pivot rule proved to be the fastest and the most robust on
780 /// various test inputs. Thus it is the default option.
782 /// \return \c true if a feasible flow can be found.
783 bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
784 return init() && start(pivot_rule);
789 /// \name Query Functions
790 /// The results of the algorithm can be obtained using these
792 /// \ref lemon::NetworkSimplex::run() "run()" must be called before
797 /// \brief Return a const reference to the edge map storing the
800 /// Return a const reference to the edge map storing the found flow.
802 /// \pre \ref run() must be called before using this function.
803 const FlowMap& flowMap() const {
804 return *_flow_result;
807 /// \brief Return a const reference to the node map storing the
808 /// found potentials (the dual solution).
810 /// Return a const reference to the node map storing the found
811 /// potentials (the dual solution).
813 /// \pre \ref run() must be called before using this function.
814 const PotentialMap& potentialMap() const {
815 return *_potential_result;
818 /// \brief Return the flow on the given edge.
820 /// Return the flow on the given edge.
822 /// \pre \ref run() must be called before using this function.
823 Capacity flow(const typename Graph::Edge& edge) const {
824 return (*_flow_result)[edge];
827 /// \brief Return the potential of the given node.
829 /// Return the potential of the given node.
831 /// \pre \ref run() must be called before using this function.
832 Cost potential(const typename Graph::Node& node) const {
833 return (*_potential_result)[node];
836 /// \brief Return the total cost of the found flow.
838 /// Return the total cost of the found flow. The complexity of the
839 /// function is \f$ O(e) \f$.
841 /// \pre \ref run() must be called before using this function.
842 Cost totalCost() const {
844 for (EdgeIt e(_orig_graph); e != INVALID; ++e)
845 c += (*_flow_result)[e] * _orig_cost[e];
853 // Initialize internal data structures
855 // Initialize result maps
857 _flow_result = new FlowMap(_orig_graph);
860 if (!_potential_result) {
861 _potential_result = new PotentialMap(_orig_graph);
862 _local_potential = true;
865 // Initialize vectors
866 _node_num = countNodes(_orig_graph);
867 _edge_num = countEdges(_orig_graph);
868 int all_node_num = _node_num + 1;
869 int all_edge_num = _edge_num + _node_num;
871 _edge.resize(_edge_num);
872 _node.reserve(_node_num);
873 _source.resize(all_edge_num);
874 _target.resize(all_edge_num);
876 _cap.resize(all_edge_num);
877 _cost.resize(all_edge_num);
878 _supply.resize(all_node_num);
879 _flow.resize(all_edge_num, 0);
880 _pi.resize(all_node_num, 0);
882 _parent.resize(all_node_num);
883 _pred.resize(all_node_num);
884 _forward.resize(all_node_num);
885 _thread.resize(all_node_num);
886 _rev_thread.resize(all_node_num);
887 _succ_num.resize(all_node_num);
888 _last_succ.resize(all_node_num);
889 _state.resize(all_edge_num, STATE_LOWER);
891 // Initialize node related data
892 bool valid_supply = true;
896 for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
899 _supply[i] = (*_orig_supply)[n];
902 valid_supply = (sum == 0);
905 for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
910 _supply[_node_id[_orig_source]] = _orig_flow_value;
911 _supply[_node_id[_orig_target]] = -_orig_flow_value;
913 if (!valid_supply) return false;
915 // Set data for the artificial root node
920 _rev_thread[0] = _root;
921 _succ_num[_root] = all_node_num;
922 _last_succ[_root] = _root - 1;
928 for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
932 // Initialize edge maps
933 for (int i = 0; i != _edge_num; ++i) {
935 _source[i] = _node_id[_orig_graph.source(e)];
936 _target[i] = _node_id[_orig_graph.target(e)];
937 _cost[i] = _orig_cost[e];
938 _cap[i] = _orig_cap[e];
941 // Remove non-zero lower bounds
943 for (int i = 0; i != _edge_num; ++i) {
944 Capacity c = (*_orig_lower)[_edge[i]];
947 _supply[_source[i]] -= c;
948 _supply[_target[i]] += c;
953 // Add artificial edges and initialize the spanning tree data structure
954 Cost max_cost = std::numeric_limits<Cost>::max() / 2 + 1;
955 Capacity max_cap = std::numeric_limits<Capacity>::max();
956 for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
960 _rev_thread[u + 1] = u;
964 _state[e] = STATE_TREE;
965 if (_supply[u] >= 0) {
970 _flow[e] = _supply[u];
978 _flow[e] = -_supply[u];
986 // Find the join node
987 void findJoinNode() {
988 int u = _source[_in_edge];
989 int v = _target[_in_edge];
991 if (_succ_num[u] < _succ_num[v]) {
1000 // Find the leaving edge of the cycle and returns true if the
1001 // leaving edge is not the same as the entering edge
1002 bool findLeavingEdge() {
1003 // Initialize first and second nodes according to the direction
1005 if (_state[_in_edge] == STATE_LOWER) {
1006 first = _source[_in_edge];
1007 second = _target[_in_edge];
1009 first = _target[_in_edge];
1010 second = _source[_in_edge];
1012 delta = _cap[_in_edge];
1017 // Search the cycle along the path form the first node to the root
1018 for (int u = first; u != join; u = _parent[u]) {
1020 d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1027 // Search the cycle along the path form the second node to the root
1028 for (int u = second; u != join; u = _parent[u]) {
1030 d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1048 // Change _flow and _state vectors
1049 void changeFlow(bool change) {
1050 // Augment along the cycle
1052 Capacity val = _state[_in_edge] * delta;
1053 _flow[_in_edge] += val;
1054 for (int u = _source[_in_edge]; u != join; u = _parent[u]) {
1055 _flow[_pred[u]] += _forward[u] ? -val : val;
1057 for (int u = _target[_in_edge]; u != join; u = _parent[u]) {
1058 _flow[_pred[u]] += _forward[u] ? val : -val;
1061 // Update the state of the entering and leaving edges
1063 _state[_in_edge] = STATE_TREE;
1064 _state[_pred[u_out]] =
1065 (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1067 _state[_in_edge] = -_state[_in_edge];
1071 // Update the tree structure
1072 void updateTreeStructure() {
1074 int old_rev_thread = _rev_thread[u_out];
1075 int old_succ_num = _succ_num[u_out];
1076 int old_last_succ = _last_succ[u_out];
1077 v_out = _parent[u_out];
1079 u = _last_succ[u_in]; // the last successor of u_in
1080 right = _thread[u]; // the node after it
1082 // Handle the case when old_rev_thread equals to v_in
1083 // (it also means that join and v_out coincide)
1084 if (old_rev_thread == v_in) {
1085 last = _thread[_last_succ[u_out]];
1087 last = _thread[v_in];
1090 // Update _thread and _parent along the stem nodes (i.e. the nodes
1091 // between u_in and u_out, whose parent have to be changed)
1092 _thread[v_in] = stem = u_in;
1093 _dirty_revs.clear();
1094 _dirty_revs.push_back(v_in);
1096 while (stem != u_out) {
1097 // Insert the next stem node into the thread list
1098 new_stem = _parent[stem];
1099 _thread[u] = new_stem;
1100 _dirty_revs.push_back(u);
1102 // Remove the subtree of stem from the thread list
1103 w = _rev_thread[stem];
1105 _rev_thread[right] = w;
1107 // Change the parent node and shift stem nodes
1108 _parent[stem] = par_stem;
1112 // Update u and right nodes
1113 u = _last_succ[stem] == _last_succ[par_stem] ?
1114 _rev_thread[par_stem] : _last_succ[stem];
1117 _parent[u_out] = par_stem;
1118 _last_succ[u_out] = u;
1120 _rev_thread[last] = u;
1122 // Remove the subtree of u_out from the thread list except for
1123 // the case when old_rev_thread equals to v_in
1124 // (it also means that join and v_out coincide)
1125 if (old_rev_thread != v_in) {
1126 _thread[old_rev_thread] = right;
1127 _rev_thread[right] = old_rev_thread;
1130 // Update _rev_thread using the new _thread values
1131 for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1133 _rev_thread[_thread[u]] = u;
1136 // Update _last_succ for the stem nodes from u_out to u_in
1137 for (u = u_out; u != u_in; u = _parent[u]) {
1138 _last_succ[_parent[u]] = _last_succ[u];
1141 // Set limits for updating _last_succ form v_in and v_out
1143 int up_limit_in = -1;
1144 int up_limit_out = -1;
1145 if (_last_succ[join] == v_in) {
1146 up_limit_out = join;
1151 // Update _last_succ from v_in towards the root
1152 for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1154 _last_succ[u] = _last_succ[u_out];
1156 // Update _last_succ from v_out towards the root
1157 if (join != old_rev_thread && v_in != old_rev_thread) {
1158 for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1160 _last_succ[u] = old_rev_thread;
1163 for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1165 _last_succ[u] = _last_succ[u_out];
1169 // Update _pred and _forward for the stem nodes from u_out to u_in
1173 _pred[u] = _pred[w];
1174 _forward[u] = !_forward[w];
1177 _pred[u_in] = _in_edge;
1178 _forward[u_in] = (u_in == _source[_in_edge]);
1180 // Update _succ_num from v_in to join
1181 for (u = v_in; u != join; u = _parent[u]) {
1182 _succ_num[u] += old_succ_num;
1184 // Update _succ_num from v_out to join
1185 for (u = v_out; u != join; u = _parent[u]) {
1186 _succ_num[u] -= old_succ_num;
1188 // Update _succ_num for the stem nodes from u_out to u_in
1193 tmp = _succ_num[u] - _succ_num[w] + tmp;
1197 _succ_num[u_in] = old_succ_num;
1200 // Update potentials
1201 void updatePotential() {
1202 Cost sigma = _forward[u_in] ?
1203 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1204 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1205 // Update in the lower subtree (which has been moved)
1206 int end = _thread[_last_succ[u_in]];
1207 for (int u = u_in; u != end; u = _thread[u]) {
1212 // Execute the algorithm
1213 bool start(PivotRuleEnum pivot_rule) {
1214 // Select the pivot rule implementation
1215 switch (pivot_rule) {
1216 case FIRST_ELIGIBLE_PIVOT:
1217 return start<FirstEligiblePivotRule>();
1218 case BEST_ELIGIBLE_PIVOT:
1219 return start<BestEligiblePivotRule>();
1220 case BLOCK_SEARCH_PIVOT:
1221 return start<BlockSearchPivotRule>();
1222 case CANDIDATE_LIST_PIVOT:
1223 return start<CandidateListPivotRule>();
1224 case ALTERING_LIST_PIVOT:
1225 return start<AlteringListPivotRule>();
1230 template<class PivotRuleImplementation>
1232 PivotRuleImplementation pivot(*this);
1234 // Execute the network simplex algorithm
1235 while (pivot.findEnteringEdge()) {
1237 bool change = findLeavingEdge();
1240 updateTreeStructure();
1245 // Check if the flow amount equals zero on all the artificial edges
1246 for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
1247 if (_flow[e] > 0) return false;
1250 // Copy flow values to _flow_result
1252 for (int i = 0; i != _edge_num; ++i) {
1254 (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
1257 for (int i = 0; i != _edge_num; ++i) {
1258 (*_flow_result)[_edge[i]] = _flow[i];
1261 // Copy potential values to _potential_result
1262 for (int i = 0; i != _node_num; ++i) {
1263 (*_potential_result)[_node[i]] = _pi[i];
1269 }; //class NetworkSimplex
1275 #endif //LEMON_NETWORK_SIMPLEX_H