47 /// |
47 /// |
48 /// In general this class is the fastest implementation available |
48 /// In general this class is the fastest implementation available |
49 /// in LEMON for the minimum cost flow problem. |
49 /// in LEMON for the minimum cost flow problem. |
50 /// |
50 /// |
51 /// \tparam GR The digraph type the algorithm runs on. |
51 /// \tparam GR The digraph type the algorithm runs on. |
52 /// \tparam V The value type used in the algorithm. |
52 /// \tparam F The value type used for flow amounts, capacity bounds |
53 /// By default it is \c int. |
53 /// and supply values in the algorithm. By default it is \c int. |
|
54 /// \tparam C The value type used for costs and potentials in the |
|
55 /// algorithm. By default it is the same as \c F. |
54 /// |
56 /// |
55 /// \warning The value type must be a signed integer type. |
57 /// \warning Both value types must be signed integer types. |
56 /// |
58 /// |
57 /// \note %NetworkSimplex provides five different pivot rule |
59 /// \note %NetworkSimplex provides five different pivot rule |
58 /// implementations. For more information see \ref PivotRule. |
60 /// implementations. For more information see \ref PivotRule. |
59 template <typename GR, typename V = int> |
61 template <typename GR, typename F = int, typename C = F> |
60 class NetworkSimplex |
62 class NetworkSimplex |
61 { |
63 { |
62 public: |
64 public: |
63 |
65 |
64 /// The value type of the algorithm |
66 /// The flow type of the algorithm |
65 typedef V Value; |
67 typedef F Flow; |
|
68 /// The cost type of the algorithm |
|
69 typedef C Cost; |
66 /// The type of the flow map |
70 /// The type of the flow map |
67 typedef typename GR::template ArcMap<Value> FlowMap; |
71 typedef typename GR::template ArcMap<Flow> FlowMap; |
68 /// The type of the potential map |
72 /// The type of the potential map |
69 typedef typename GR::template NodeMap<Value> PotentialMap; |
73 typedef typename GR::template NodeMap<Cost> PotentialMap; |
70 |
74 |
71 public: |
75 public: |
72 |
76 |
73 /// \brief Enum type for selecting the pivot rule. |
77 /// \brief Enum type for selecting the pivot rule. |
74 /// |
78 /// |
115 |
119 |
116 private: |
120 private: |
117 |
121 |
118 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
122 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
119 |
123 |
120 typedef typename GR::template ArcMap<Value> ValueArcMap; |
124 typedef typename GR::template ArcMap<Flow> FlowArcMap; |
121 typedef typename GR::template NodeMap<Value> ValueNodeMap; |
125 typedef typename GR::template ArcMap<Cost> CostArcMap; |
|
126 typedef typename GR::template NodeMap<Flow> FlowNodeMap; |
122 |
127 |
123 typedef std::vector<Arc> ArcVector; |
128 typedef std::vector<Arc> ArcVector; |
124 typedef std::vector<Node> NodeVector; |
129 typedef std::vector<Node> NodeVector; |
125 typedef std::vector<int> IntVector; |
130 typedef std::vector<int> IntVector; |
126 typedef std::vector<bool> BoolVector; |
131 typedef std::vector<bool> BoolVector; |
127 typedef std::vector<Value> ValueVector; |
132 typedef std::vector<Flow> FlowVector; |
|
133 typedef std::vector<Cost> CostVector; |
128 |
134 |
129 // State constants for arcs |
135 // State constants for arcs |
130 enum ArcStateEnum { |
136 enum ArcStateEnum { |
131 STATE_UPPER = -1, |
137 STATE_UPPER = -1, |
132 STATE_TREE = 0, |
138 STATE_TREE = 0, |
139 const GR &_graph; |
145 const GR &_graph; |
140 int _node_num; |
146 int _node_num; |
141 int _arc_num; |
147 int _arc_num; |
142 |
148 |
143 // Parameters of the problem |
149 // Parameters of the problem |
144 ValueArcMap *_plower; |
150 FlowArcMap *_plower; |
145 ValueArcMap *_pupper; |
151 FlowArcMap *_pupper; |
146 ValueArcMap *_pcost; |
152 CostArcMap *_pcost; |
147 ValueNodeMap *_psupply; |
153 FlowNodeMap *_psupply; |
148 bool _pstsup; |
154 bool _pstsup; |
149 Node _psource, _ptarget; |
155 Node _psource, _ptarget; |
150 Value _pstflow; |
156 Flow _pstflow; |
151 |
157 |
152 // Result maps |
158 // Result maps |
153 FlowMap *_flow_map; |
159 FlowMap *_flow_map; |
154 PotentialMap *_potential_map; |
160 PotentialMap *_potential_map; |
155 bool _local_flow; |
161 bool _local_flow; |
182 |
188 |
183 // Temporary data used in the current pivot iteration |
189 // Temporary data used in the current pivot iteration |
184 int in_arc, join, u_in, v_in, u_out, v_out; |
190 int in_arc, join, u_in, v_in, u_out, v_out; |
185 int first, second, right, last; |
191 int first, second, right, last; |
186 int stem, par_stem, new_stem; |
192 int stem, par_stem, new_stem; |
187 Value delta; |
193 Flow delta; |
188 |
194 |
189 private: |
195 private: |
190 |
196 |
191 // Implementation of the First Eligible pivot rule |
197 // Implementation of the First Eligible pivot rule |
192 class FirstEligiblePivotRule |
198 class FirstEligiblePivotRule |
194 private: |
200 private: |
195 |
201 |
196 // References to the NetworkSimplex class |
202 // References to the NetworkSimplex class |
197 const IntVector &_source; |
203 const IntVector &_source; |
198 const IntVector &_target; |
204 const IntVector &_target; |
199 const ValueVector &_cost; |
205 const CostVector &_cost; |
200 const IntVector &_state; |
206 const IntVector &_state; |
201 const ValueVector &_pi; |
207 const CostVector &_pi; |
202 int &_in_arc; |
208 int &_in_arc; |
203 int _arc_num; |
209 int _arc_num; |
204 |
210 |
205 // Pivot rule data |
211 // Pivot rule data |
206 int _next_arc; |
212 int _next_arc; |
214 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
220 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
215 {} |
221 {} |
216 |
222 |
217 // Find next entering arc |
223 // Find next entering arc |
218 bool findEnteringArc() { |
224 bool findEnteringArc() { |
219 Value c; |
225 Cost c; |
220 for (int e = _next_arc; e < _arc_num; ++e) { |
226 for (int e = _next_arc; e < _arc_num; ++e) { |
221 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
227 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
222 if (c < 0) { |
228 if (c < 0) { |
223 _in_arc = e; |
229 _in_arc = e; |
224 _next_arc = e + 1; |
230 _next_arc = e + 1; |
245 private: |
251 private: |
246 |
252 |
247 // References to the NetworkSimplex class |
253 // References to the NetworkSimplex class |
248 const IntVector &_source; |
254 const IntVector &_source; |
249 const IntVector &_target; |
255 const IntVector &_target; |
250 const ValueVector &_cost; |
256 const CostVector &_cost; |
251 const IntVector &_state; |
257 const IntVector &_state; |
252 const ValueVector &_pi; |
258 const CostVector &_pi; |
253 int &_in_arc; |
259 int &_in_arc; |
254 int _arc_num; |
260 int _arc_num; |
255 |
261 |
256 public: |
262 public: |
257 |
263 |
262 _in_arc(ns.in_arc), _arc_num(ns._arc_num) |
268 _in_arc(ns.in_arc), _arc_num(ns._arc_num) |
263 {} |
269 {} |
264 |
270 |
265 // Find next entering arc |
271 // Find next entering arc |
266 bool findEnteringArc() { |
272 bool findEnteringArc() { |
267 Value c, min = 0; |
273 Cost c, min = 0; |
268 for (int e = 0; e < _arc_num; ++e) { |
274 for (int e = 0; e < _arc_num; ++e) { |
269 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
275 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
270 if (c < min) { |
276 if (c < min) { |
271 min = c; |
277 min = c; |
272 _in_arc = e; |
278 _in_arc = e; |
284 private: |
290 private: |
285 |
291 |
286 // References to the NetworkSimplex class |
292 // References to the NetworkSimplex class |
287 const IntVector &_source; |
293 const IntVector &_source; |
288 const IntVector &_target; |
294 const IntVector &_target; |
289 const ValueVector &_cost; |
295 const CostVector &_cost; |
290 const IntVector &_state; |
296 const IntVector &_state; |
291 const ValueVector &_pi; |
297 const CostVector &_pi; |
292 int &_in_arc; |
298 int &_in_arc; |
293 int _arc_num; |
299 int _arc_num; |
294 |
300 |
295 // Pivot rule data |
301 // Pivot rule data |
296 int _block_size; |
302 int _block_size; |
312 MIN_BLOCK_SIZE ); |
318 MIN_BLOCK_SIZE ); |
313 } |
319 } |
314 |
320 |
315 // Find next entering arc |
321 // Find next entering arc |
316 bool findEnteringArc() { |
322 bool findEnteringArc() { |
317 Value c, min = 0; |
323 Cost c, min = 0; |
318 int cnt = _block_size; |
324 int cnt = _block_size; |
319 int e, min_arc = _next_arc; |
325 int e, min_arc = _next_arc; |
320 for (e = _next_arc; e < _arc_num; ++e) { |
326 for (e = _next_arc; e < _arc_num; ++e) { |
321 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
327 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
322 if (c < min) { |
328 if (c < min) { |
356 private: |
362 private: |
357 |
363 |
358 // References to the NetworkSimplex class |
364 // References to the NetworkSimplex class |
359 const IntVector &_source; |
365 const IntVector &_source; |
360 const IntVector &_target; |
366 const IntVector &_target; |
361 const ValueVector &_cost; |
367 const CostVector &_cost; |
362 const IntVector &_state; |
368 const IntVector &_state; |
363 const ValueVector &_pi; |
369 const CostVector &_pi; |
364 int &_in_arc; |
370 int &_in_arc; |
365 int _arc_num; |
371 int _arc_num; |
366 |
372 |
367 // Pivot rule data |
373 // Pivot rule data |
368 IntVector _candidates; |
374 IntVector _candidates; |
392 _candidates.resize(_list_length); |
398 _candidates.resize(_list_length); |
393 } |
399 } |
394 |
400 |
395 /// Find next entering arc |
401 /// Find next entering arc |
396 bool findEnteringArc() { |
402 bool findEnteringArc() { |
397 Value min, c; |
403 Cost min, c; |
398 int e, min_arc = _next_arc; |
404 int e, min_arc = _next_arc; |
399 if (_curr_length > 0 && _minor_count < _minor_limit) { |
405 if (_curr_length > 0 && _minor_count < _minor_limit) { |
400 // Minor iteration: select the best eligible arc from the |
406 // Minor iteration: select the best eligible arc from the |
401 // current candidate list |
407 // current candidate list |
402 ++_minor_count; |
408 ++_minor_count; |
461 private: |
467 private: |
462 |
468 |
463 // References to the NetworkSimplex class |
469 // References to the NetworkSimplex class |
464 const IntVector &_source; |
470 const IntVector &_source; |
465 const IntVector &_target; |
471 const IntVector &_target; |
466 const ValueVector &_cost; |
472 const CostVector &_cost; |
467 const IntVector &_state; |
473 const IntVector &_state; |
468 const ValueVector &_pi; |
474 const CostVector &_pi; |
469 int &_in_arc; |
475 int &_in_arc; |
470 int _arc_num; |
476 int _arc_num; |
471 |
477 |
472 // Pivot rule data |
478 // Pivot rule data |
473 int _block_size, _head_length, _curr_length; |
479 int _block_size, _head_length, _curr_length; |
474 int _next_arc; |
480 int _next_arc; |
475 IntVector _candidates; |
481 IntVector _candidates; |
476 ValueVector _cand_cost; |
482 CostVector _cand_cost; |
477 |
483 |
478 // Functor class to compare arcs during sort of the candidate list |
484 // Functor class to compare arcs during sort of the candidate list |
479 class SortFunc |
485 class SortFunc |
480 { |
486 { |
481 private: |
487 private: |
482 const ValueVector &_map; |
488 const CostVector &_map; |
483 public: |
489 public: |
484 SortFunc(const ValueVector &map) : _map(map) {} |
490 SortFunc(const CostVector &map) : _map(map) {} |
485 bool operator()(int left, int right) { |
491 bool operator()(int left, int right) { |
486 return _map[left] > _map[right]; |
492 return _map[left] > _map[right]; |
487 } |
493 } |
488 }; |
494 }; |
489 |
495 |
588 _psupply(NULL), _pstsup(false), |
594 _psupply(NULL), _pstsup(false), |
589 _flow_map(NULL), _potential_map(NULL), |
595 _flow_map(NULL), _potential_map(NULL), |
590 _local_flow(false), _local_potential(false), |
596 _local_flow(false), _local_potential(false), |
591 _node_id(graph) |
597 _node_id(graph) |
592 { |
598 { |
593 LEMON_ASSERT(std::numeric_limits<Value>::is_integer && |
599 LEMON_ASSERT(std::numeric_limits<Flow>::is_integer && |
594 std::numeric_limits<Value>::is_signed, |
600 std::numeric_limits<Flow>::is_signed, |
595 "The value type of NetworkSimplex must be a signed integer"); |
601 "The flow type of NetworkSimplex must be signed integer"); |
|
602 LEMON_ASSERT(std::numeric_limits<Cost>::is_integer && |
|
603 std::numeric_limits<Cost>::is_signed, |
|
604 "The cost type of NetworkSimplex must be signed integer"); |
596 } |
605 } |
597 |
606 |
598 /// Destructor. |
607 /// Destructor. |
599 ~NetworkSimplex() { |
608 ~NetworkSimplex() { |
600 if (_local_flow) delete _flow_map; |
609 if (_local_flow) delete _flow_map; |
607 /// If neither this function nor \ref boundMaps() is used before |
616 /// If neither this function nor \ref boundMaps() is used before |
608 /// calling \ref run(), the lower bounds will be set to zero |
617 /// calling \ref run(), the lower bounds will be set to zero |
609 /// on all arcs. |
618 /// on all arcs. |
610 /// |
619 /// |
611 /// \param map An arc map storing the lower bounds. |
620 /// \param map An arc map storing the lower bounds. |
612 /// Its \c Value type must be convertible to the \c Value type |
621 /// Its \c Value type must be convertible to the \c Flow type |
613 /// of the algorithm. |
622 /// of the algorithm. |
614 /// |
623 /// |
615 /// \return <tt>(*this)</tt> |
624 /// \return <tt>(*this)</tt> |
616 template <typename LOWER> |
625 template <typename LOWER> |
617 NetworkSimplex& lowerMap(const LOWER& map) { |
626 NetworkSimplex& lowerMap(const LOWER& map) { |
618 delete _plower; |
627 delete _plower; |
619 _plower = new ValueArcMap(_graph); |
628 _plower = new FlowArcMap(_graph); |
620 for (ArcIt a(_graph); a != INVALID; ++a) { |
629 for (ArcIt a(_graph); a != INVALID; ++a) { |
621 (*_plower)[a] = map[a]; |
630 (*_plower)[a] = map[a]; |
622 } |
631 } |
623 return *this; |
632 return *this; |
624 } |
633 } |
627 /// |
636 /// |
628 /// This function sets the upper bounds (capacities) on the arcs. |
637 /// This function sets the upper bounds (capacities) on the arcs. |
629 /// If none of the functions \ref upperMap(), \ref capacityMap() |
638 /// If none of the functions \ref upperMap(), \ref capacityMap() |
630 /// and \ref boundMaps() is used before calling \ref run(), |
639 /// and \ref boundMaps() is used before calling \ref run(), |
631 /// the upper bounds (capacities) will be set to |
640 /// the upper bounds (capacities) will be set to |
632 /// \c std::numeric_limits<Value>::max() on all arcs. |
641 /// \c std::numeric_limits<Flow>::max() on all arcs. |
633 /// |
642 /// |
634 /// \param map An arc map storing the upper bounds. |
643 /// \param map An arc map storing the upper bounds. |
635 /// Its \c Value type must be convertible to the \c Value type |
644 /// Its \c Value type must be convertible to the \c Flow type |
636 /// of the algorithm. |
645 /// of the algorithm. |
637 /// |
646 /// |
638 /// \return <tt>(*this)</tt> |
647 /// \return <tt>(*this)</tt> |
639 template<typename UPPER> |
648 template<typename UPPER> |
640 NetworkSimplex& upperMap(const UPPER& map) { |
649 NetworkSimplex& upperMap(const UPPER& map) { |
641 delete _pupper; |
650 delete _pupper; |
642 _pupper = new ValueArcMap(_graph); |
651 _pupper = new FlowArcMap(_graph); |
643 for (ArcIt a(_graph); a != INVALID; ++a) { |
652 for (ArcIt a(_graph); a != INVALID; ++a) { |
644 (*_pupper)[a] = map[a]; |
653 (*_pupper)[a] = map[a]; |
645 } |
654 } |
646 return *this; |
655 return *this; |
647 } |
656 } |
664 /// calling \ref run(), the lower bounds will be set to zero |
673 /// calling \ref run(), the lower bounds will be set to zero |
665 /// on all arcs. |
674 /// on all arcs. |
666 /// If none of the functions \ref upperMap(), \ref capacityMap() |
675 /// If none of the functions \ref upperMap(), \ref capacityMap() |
667 /// and \ref boundMaps() is used before calling \ref run(), |
676 /// and \ref boundMaps() is used before calling \ref run(), |
668 /// the upper bounds (capacities) will be set to |
677 /// the upper bounds (capacities) will be set to |
669 /// \c std::numeric_limits<Value>::max() on all arcs. |
678 /// \c std::numeric_limits<Flow>::max() on all arcs. |
670 /// |
679 /// |
671 /// \param lower An arc map storing the lower bounds. |
680 /// \param lower An arc map storing the lower bounds. |
672 /// \param upper An arc map storing the upper bounds. |
681 /// \param upper An arc map storing the upper bounds. |
673 /// |
682 /// |
674 /// The \c Value type of the maps must be convertible to the |
683 /// The \c Value type of the maps must be convertible to the |
675 /// \c Value type of the algorithm. |
684 /// \c Flow type of the algorithm. |
676 /// |
685 /// |
677 /// \note This function is just a shortcut of calling \ref lowerMap() |
686 /// \note This function is just a shortcut of calling \ref lowerMap() |
678 /// and \ref upperMap() separately. |
687 /// and \ref upperMap() separately. |
679 /// |
688 /// |
680 /// \return <tt>(*this)</tt> |
689 /// \return <tt>(*this)</tt> |
688 /// This function sets the costs of the arcs. |
697 /// This function sets the costs of the arcs. |
689 /// If it is not used before calling \ref run(), the costs |
698 /// If it is not used before calling \ref run(), the costs |
690 /// will be set to \c 1 on all arcs. |
699 /// will be set to \c 1 on all arcs. |
691 /// |
700 /// |
692 /// \param map An arc map storing the costs. |
701 /// \param map An arc map storing the costs. |
693 /// Its \c Value type must be convertible to the \c Value type |
702 /// Its \c Value type must be convertible to the \c Cost type |
694 /// of the algorithm. |
703 /// of the algorithm. |
695 /// |
704 /// |
696 /// \return <tt>(*this)</tt> |
705 /// \return <tt>(*this)</tt> |
697 template<typename COST> |
706 template<typename COST> |
698 NetworkSimplex& costMap(const COST& map) { |
707 NetworkSimplex& costMap(const COST& map) { |
699 delete _pcost; |
708 delete _pcost; |
700 _pcost = new ValueArcMap(_graph); |
709 _pcost = new CostArcMap(_graph); |
701 for (ArcIt a(_graph); a != INVALID; ++a) { |
710 for (ArcIt a(_graph); a != INVALID; ++a) { |
702 (*_pcost)[a] = map[a]; |
711 (*_pcost)[a] = map[a]; |
703 } |
712 } |
704 return *this; |
713 return *this; |
705 } |
714 } |
710 /// If neither this function nor \ref stSupply() is used before |
719 /// If neither this function nor \ref stSupply() is used before |
711 /// calling \ref run(), the supply of each node will be set to zero. |
720 /// calling \ref run(), the supply of each node will be set to zero. |
712 /// (It makes sense only if non-zero lower bounds are given.) |
721 /// (It makes sense only if non-zero lower bounds are given.) |
713 /// |
722 /// |
714 /// \param map A node map storing the supply values. |
723 /// \param map A node map storing the supply values. |
715 /// Its \c Value type must be convertible to the \c Value type |
724 /// Its \c Value type must be convertible to the \c Flow type |
716 /// of the algorithm. |
725 /// of the algorithm. |
717 /// |
726 /// |
718 /// \return <tt>(*this)</tt> |
727 /// \return <tt>(*this)</tt> |
719 template<typename SUP> |
728 template<typename SUP> |
720 NetworkSimplex& supplyMap(const SUP& map) { |
729 NetworkSimplex& supplyMap(const SUP& map) { |
721 delete _psupply; |
730 delete _psupply; |
722 _pstsup = false; |
731 _pstsup = false; |
723 _psupply = new ValueNodeMap(_graph); |
732 _psupply = new FlowNodeMap(_graph); |
724 for (NodeIt n(_graph); n != INVALID; ++n) { |
733 for (NodeIt n(_graph); n != INVALID; ++n) { |
725 (*_psupply)[n] = map[n]; |
734 (*_psupply)[n] = map[n]; |
726 } |
735 } |
727 return *this; |
736 return *this; |
728 } |
737 } |
739 /// \param t The target node. |
748 /// \param t The target node. |
740 /// \param k The required amount of flow from node \c s to node \c t |
749 /// \param k The required amount of flow from node \c s to node \c t |
741 /// (i.e. the supply of \c s and the demand of \c t). |
750 /// (i.e. the supply of \c s and the demand of \c t). |
742 /// |
751 /// |
743 /// \return <tt>(*this)</tt> |
752 /// \return <tt>(*this)</tt> |
744 NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { |
753 NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) { |
745 delete _psupply; |
754 delete _psupply; |
746 _psupply = NULL; |
755 _psupply = NULL; |
747 _pstsup = true; |
756 _pstsup = true; |
748 _psource = s; |
757 _psource = s; |
749 _ptarget = t; |
758 _ptarget = t; |
872 /// @{ |
881 /// @{ |
873 |
882 |
874 /// \brief Return the total cost of the found flow. |
883 /// \brief Return the total cost of the found flow. |
875 /// |
884 /// |
876 /// This function returns the total cost of the found flow. |
885 /// This function returns the total cost of the found flow. |
877 /// The complexity of the function is \f$ O(e) \f$. |
886 /// The complexity of the function is O(e). |
878 /// |
887 /// |
879 /// \note The return type of the function can be specified as a |
888 /// \note The return type of the function can be specified as a |
880 /// template parameter. For example, |
889 /// template parameter. For example, |
881 /// \code |
890 /// \code |
882 /// ns.totalCost<double>(); |
891 /// ns.totalCost<double>(); |
883 /// \endcode |
892 /// \endcode |
884 /// It is useful if the total cost cannot be stored in the \c Value |
893 /// It is useful if the total cost cannot be stored in the \c Cost |
885 /// type of the algorithm, which is the default return type of the |
894 /// type of the algorithm, which is the default return type of the |
886 /// function. |
895 /// function. |
887 /// |
896 /// |
888 /// \pre \ref run() must be called before using this function. |
897 /// \pre \ref run() must be called before using this function. |
889 template <typename Num> |
898 template <typename Num> |
898 } |
907 } |
899 return c; |
908 return c; |
900 } |
909 } |
901 |
910 |
902 #ifndef DOXYGEN |
911 #ifndef DOXYGEN |
903 Value totalCost() const { |
912 Cost totalCost() const { |
904 return totalCost<Value>(); |
913 return totalCost<Cost>(); |
905 } |
914 } |
906 #endif |
915 #endif |
907 |
916 |
908 /// \brief Return the flow on the given arc. |
917 /// \brief Return the flow on the given arc. |
909 /// |
918 /// |
910 /// This function returns the flow on the given arc. |
919 /// This function returns the flow on the given arc. |
911 /// |
920 /// |
912 /// \pre \ref run() must be called before using this function. |
921 /// \pre \ref run() must be called before using this function. |
913 Value flow(const Arc& a) const { |
922 Flow flow(const Arc& a) const { |
914 return (*_flow_map)[a]; |
923 return (*_flow_map)[a]; |
915 } |
924 } |
916 |
925 |
917 /// \brief Return a const reference to the flow map. |
926 /// \brief Return a const reference to the flow map. |
918 /// |
927 /// |
928 /// |
937 /// |
929 /// This function returns the potential (dual value) of the |
938 /// This function returns the potential (dual value) of the |
930 /// given node. |
939 /// given node. |
931 /// |
940 /// |
932 /// \pre \ref run() must be called before using this function. |
941 /// \pre \ref run() must be called before using this function. |
933 Value potential(const Node& n) const { |
942 Cost potential(const Node& n) const { |
934 return (*_potential_map)[n]; |
943 return (*_potential_map)[n]; |
935 } |
944 } |
936 |
945 |
937 /// \brief Return a const reference to the potential map |
946 /// \brief Return a const reference to the potential map |
938 /// (the dual solution). |
947 /// (the dual solution). |
1033 _arc_ref[i] = e; |
1042 _arc_ref[i] = e; |
1034 if ((i += k) >= _arc_num) i = (i % k) + 1; |
1043 if ((i += k) >= _arc_num) i = (i % k) + 1; |
1035 } |
1044 } |
1036 |
1045 |
1037 // Initialize arc maps |
1046 // Initialize arc maps |
|
1047 Flow max_cap = std::numeric_limits<Flow>::max(); |
|
1048 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
1038 if (_pupper && _pcost) { |
1049 if (_pupper && _pcost) { |
1039 for (int i = 0; i != _arc_num; ++i) { |
1050 for (int i = 0; i != _arc_num; ++i) { |
1040 Arc e = _arc_ref[i]; |
1051 Arc e = _arc_ref[i]; |
1041 _source[i] = _node_id[_graph.source(e)]; |
1052 _source[i] = _node_id[_graph.source(e)]; |
1042 _target[i] = _node_id[_graph.target(e)]; |
1053 _target[i] = _node_id[_graph.target(e)]; |
1055 } |
1066 } |
1056 if (_pupper) { |
1067 if (_pupper) { |
1057 for (int i = 0; i != _arc_num; ++i) |
1068 for (int i = 0; i != _arc_num; ++i) |
1058 _cap[i] = (*_pupper)[_arc_ref[i]]; |
1069 _cap[i] = (*_pupper)[_arc_ref[i]]; |
1059 } else { |
1070 } else { |
1060 Value val = std::numeric_limits<Value>::max(); |
|
1061 for (int i = 0; i != _arc_num; ++i) |
1071 for (int i = 0; i != _arc_num; ++i) |
1062 _cap[i] = val; |
1072 _cap[i] = max_cap; |
1063 } |
1073 } |
1064 if (_pcost) { |
1074 if (_pcost) { |
1065 for (int i = 0; i != _arc_num; ++i) |
1075 for (int i = 0; i != _arc_num; ++i) |
1066 _cost[i] = (*_pcost)[_arc_ref[i]]; |
1076 _cost[i] = (*_pcost)[_arc_ref[i]]; |
1067 } else { |
1077 } else { |
1071 } |
1081 } |
1072 |
1082 |
1073 // Remove non-zero lower bounds |
1083 // Remove non-zero lower bounds |
1074 if (_plower) { |
1084 if (_plower) { |
1075 for (int i = 0; i != _arc_num; ++i) { |
1085 for (int i = 0; i != _arc_num; ++i) { |
1076 Value c = (*_plower)[_arc_ref[i]]; |
1086 Flow c = (*_plower)[_arc_ref[i]]; |
1077 if (c != 0) { |
1087 if (c != 0) { |
1078 _cap[i] -= c; |
1088 _cap[i] -= c; |
1079 _supply[_source[i]] -= c; |
1089 _supply[_source[i]] -= c; |
1080 _supply[_target[i]] += c; |
1090 _supply[_target[i]] += c; |
1081 } |
1091 } |
1082 } |
1092 } |
1083 } |
1093 } |
1084 |
1094 |
1085 // Add artificial arcs and initialize the spanning tree data structure |
1095 // Add artificial arcs and initialize the spanning tree data structure |
1086 Value max_cap = std::numeric_limits<Value>::max(); |
|
1087 Value max_cost = std::numeric_limits<Value>::max() / 4; |
|
1088 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
1096 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
1089 _thread[u] = u + 1; |
1097 _thread[u] = u + 1; |
1090 _rev_thread[u + 1] = u; |
1098 _rev_thread[u + 1] = u; |
1091 _succ_num[u] = 1; |
1099 _succ_num[u] = 1; |
1092 _last_succ[u] = u; |
1100 _last_succ[u] = u; |
1135 first = _target[in_arc]; |
1143 first = _target[in_arc]; |
1136 second = _source[in_arc]; |
1144 second = _source[in_arc]; |
1137 } |
1145 } |
1138 delta = _cap[in_arc]; |
1146 delta = _cap[in_arc]; |
1139 int result = 0; |
1147 int result = 0; |
1140 Value d; |
1148 Flow d; |
1141 int e; |
1149 int e; |
1142 |
1150 |
1143 // Search the cycle along the path form the first node to the root |
1151 // Search the cycle along the path form the first node to the root |
1144 for (int u = first; u != join; u = _parent[u]) { |
1152 for (int u = first; u != join; u = _parent[u]) { |
1145 e = _pred[u]; |
1153 e = _pred[u]; |
1173 |
1181 |
1174 // Change _flow and _state vectors |
1182 // Change _flow and _state vectors |
1175 void changeFlow(bool change) { |
1183 void changeFlow(bool change) { |
1176 // Augment along the cycle |
1184 // Augment along the cycle |
1177 if (delta > 0) { |
1185 if (delta > 0) { |
1178 Value val = _state[in_arc] * delta; |
1186 Flow val = _state[in_arc] * delta; |
1179 _flow[in_arc] += val; |
1187 _flow[in_arc] += val; |
1180 for (int u = _source[in_arc]; u != join; u = _parent[u]) { |
1188 for (int u = _source[in_arc]; u != join; u = _parent[u]) { |
1181 _flow[_pred[u]] += _forward[u] ? -val : val; |
1189 _flow[_pred[u]] += _forward[u] ? -val : val; |
1182 } |
1190 } |
1183 for (int u = _target[in_arc]; u != join; u = _parent[u]) { |
1191 for (int u = _target[in_arc]; u != join; u = _parent[u]) { |
1314 } |
1322 } |
1315 } |
1323 } |
1316 |
1324 |
1317 // Update potentials |
1325 // Update potentials |
1318 void updatePotential() { |
1326 void updatePotential() { |
1319 Value sigma = _forward[u_in] ? |
1327 Cost sigma = _forward[u_in] ? |
1320 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1328 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1321 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1329 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1322 if (_succ_num[u_in] > _node_num / 2) { |
1330 if (_succ_num[u_in] > _node_num / 2) { |
1323 // Update in the upper subtree (which contains the root) |
1331 // Update in the upper subtree (which contains the root) |
1324 int before = _rev_thread[u_in]; |
1332 int before = _rev_thread[u_in]; |