changeset 605 | 5232721b3f14 |
parent 604 | 8c3112a66878 |
child 606 | c7d160f73d52 |
2:5f9d9ef3182e | 3:6a4d274c4740 |
---|---|
20 #define LEMON_NETWORK_SIMPLEX_H |
20 #define LEMON_NETWORK_SIMPLEX_H |
21 |
21 |
22 /// \ingroup min_cost_flow |
22 /// \ingroup min_cost_flow |
23 /// |
23 /// |
24 /// \file |
24 /// \file |
25 /// \brief Network simplex algorithm for finding a minimum cost flow. |
25 /// \brief Network Simplex algorithm for finding a minimum cost flow. |
26 |
26 |
27 #include <vector> |
27 #include <vector> |
28 #include <limits> |
28 #include <limits> |
29 #include <algorithm> |
29 #include <algorithm> |
30 |
30 |
34 namespace lemon { |
34 namespace lemon { |
35 |
35 |
36 /// \addtogroup min_cost_flow |
36 /// \addtogroup min_cost_flow |
37 /// @{ |
37 /// @{ |
38 |
38 |
39 /// \brief Implementation of the primal network simplex algorithm |
39 /// \brief Implementation of the primal Network Simplex algorithm |
40 /// for finding a \ref min_cost_flow "minimum cost flow". |
40 /// for finding a \ref min_cost_flow "minimum cost flow". |
41 /// |
41 /// |
42 /// \ref NetworkSimplex implements the primal network simplex algorithm |
42 /// \ref NetworkSimplex implements the primal Network Simplex algorithm |
43 /// for finding a \ref min_cost_flow "minimum cost flow". |
43 /// for finding a \ref min_cost_flow "minimum cost flow". |
44 /// |
44 /// |
45 /// \tparam Digraph The digraph type the algorithm runs on. |
45 /// \tparam GR The digraph type the algorithm runs on. |
46 /// \tparam LowerMap The type of the lower bound map. |
46 /// \tparam V The value type used in the algorithm. |
47 /// \tparam CapacityMap The type of the capacity (upper bound) map. |
47 /// By default it is \c int. |
48 /// \tparam CostMap The type of the cost (length) map. |
|
49 /// \tparam SupplyMap The type of the supply map. |
|
50 /// |
48 /// |
51 /// \warning |
49 /// \warning \c V must be a signed integer type. |
52 /// - Arc 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. |
|
56 /// |
50 /// |
57 /// \note \ref NetworkSimplex provides five different pivot rule |
51 /// \note %NetworkSimplex provides five different pivot rule |
58 /// implementations that significantly affect the efficiency of the |
52 /// implementations. For more information see \ref PivotRule. |
59 /// algorithm. |
53 template <typename GR, typename V = int> |
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. |
|
64 #ifdef DOXYGEN |
|
65 template < typename Digraph, |
|
66 typename LowerMap, |
|
67 typename CapacityMap, |
|
68 typename CostMap, |
|
69 typename SupplyMap > |
|
70 |
|
71 #else |
|
72 template < typename Digraph, |
|
73 typename LowerMap = typename Digraph::template ArcMap<int>, |
|
74 typename CapacityMap = typename Digraph::template ArcMap<int>, |
|
75 typename CostMap = typename Digraph::template ArcMap<int>, |
|
76 typename SupplyMap = typename Digraph::template NodeMap<int> > |
|
77 #endif |
|
78 class NetworkSimplex |
54 class NetworkSimplex |
79 { |
55 { |
80 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
56 public: |
81 |
57 |
82 typedef typename CapacityMap::Value Capacity; |
58 /// The value type of the algorithm |
83 typedef typename CostMap::Value Cost; |
59 typedef V Value; |
84 typedef typename SupplyMap::Value Supply; |
60 /// The type of the flow map |
61 typedef typename GR::template ArcMap<Value> FlowMap; |
|
62 /// The type of the potential map |
|
63 typedef typename GR::template NodeMap<Value> PotentialMap; |
|
64 |
|
65 public: |
|
66 |
|
67 /// \brief Enum type for selecting the pivot rule. |
|
68 /// |
|
69 /// Enum type for selecting the pivot rule for the \ref run() |
|
70 /// function. |
|
71 /// |
|
72 /// \ref NetworkSimplex provides five different pivot rule |
|
73 /// implementations that significantly affect the running time |
|
74 /// of the algorithm. |
|
75 /// By default \ref BLOCK_SEARCH "Block Search" is used, which |
|
76 /// proved to be the most efficient and the most robust on various |
|
77 /// test inputs according to our benchmark tests. |
|
78 /// However another pivot rule can be selected using the \ref run() |
|
79 /// function with the proper parameter. |
|
80 enum PivotRule { |
|
81 |
|
82 /// The First Eligible pivot rule. |
|
83 /// The next eligible arc is selected in a wraparound fashion |
|
84 /// in every iteration. |
|
85 FIRST_ELIGIBLE, |
|
86 |
|
87 /// The Best Eligible pivot rule. |
|
88 /// The best eligible arc is selected in every iteration. |
|
89 BEST_ELIGIBLE, |
|
90 |
|
91 /// The Block Search pivot rule. |
|
92 /// A specified number of arcs are examined in every iteration |
|
93 /// in a wraparound fashion and the best eligible arc is selected |
|
94 /// from this block. |
|
95 BLOCK_SEARCH, |
|
96 |
|
97 /// The Candidate List pivot rule. |
|
98 /// In a major iteration a candidate list is built from eligible arcs |
|
99 /// in a wraparound fashion and in the following minor iterations |
|
100 /// the best eligible arc is selected from this list. |
|
101 CANDIDATE_LIST, |
|
102 |
|
103 /// The Altering Candidate List pivot rule. |
|
104 /// It is a modified version of the Candidate List method. |
|
105 /// It keeps only the several best eligible arcs from the former |
|
106 /// candidate list and extends this list in every iteration. |
|
107 ALTERING_LIST |
|
108 }; |
|
109 |
|
110 private: |
|
111 |
|
112 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
|
113 |
|
114 typedef typename GR::template ArcMap<Value> ValueArcMap; |
|
115 typedef typename GR::template NodeMap<Value> ValueNodeMap; |
|
85 |
116 |
86 typedef std::vector<Arc> ArcVector; |
117 typedef std::vector<Arc> ArcVector; |
87 typedef std::vector<Node> NodeVector; |
118 typedef std::vector<Node> NodeVector; |
88 typedef std::vector<int> IntVector; |
119 typedef std::vector<int> IntVector; |
89 typedef std::vector<bool> BoolVector; |
120 typedef std::vector<bool> BoolVector; |
90 typedef std::vector<Capacity> CapacityVector; |
121 typedef std::vector<Value> ValueVector; |
91 typedef std::vector<Cost> CostVector; |
|
92 typedef std::vector<Supply> SupplyVector; |
|
93 |
|
94 public: |
|
95 |
|
96 /// The type of the flow map |
|
97 typedef typename Digraph::template ArcMap<Capacity> FlowMap; |
|
98 /// The type of the potential map |
|
99 typedef typename Digraph::template NodeMap<Cost> PotentialMap; |
|
100 |
|
101 public: |
|
102 |
|
103 /// Enum type for selecting the pivot rule used by \ref run() |
|
104 enum PivotRuleEnum { |
|
105 FIRST_ELIGIBLE_PIVOT, |
|
106 BEST_ELIGIBLE_PIVOT, |
|
107 BLOCK_SEARCH_PIVOT, |
|
108 CANDIDATE_LIST_PIVOT, |
|
109 ALTERING_LIST_PIVOT |
|
110 }; |
|
111 |
|
112 private: |
|
113 |
122 |
114 // State constants for arcs |
123 // State constants for arcs |
115 enum ArcStateEnum { |
124 enum ArcStateEnum { |
116 STATE_UPPER = -1, |
125 STATE_UPPER = -1, |
117 STATE_TREE = 0, |
126 STATE_TREE = 0, |
118 STATE_LOWER = 1 |
127 STATE_LOWER = 1 |
119 }; |
128 }; |
120 |
129 |
121 private: |
130 private: |
122 |
131 |
123 // References for the original data |
132 // Data related to the underlying digraph |
124 const Digraph &_graph; |
133 const GR &_graph; |
125 const LowerMap *_orig_lower; |
134 int _node_num; |
126 const CapacityMap &_orig_cap; |
135 int _arc_num; |
127 const CostMap &_orig_cost; |
136 |
128 const SupplyMap *_orig_supply; |
137 // Parameters of the problem |
129 Node _orig_source; |
138 ValueArcMap *_plower; |
130 Node _orig_target; |
139 ValueArcMap *_pupper; |
131 Capacity _orig_flow_value; |
140 ValueArcMap *_pcost; |
141 ValueNodeMap *_psupply; |
|
142 bool _pstsup; |
|
143 Node _psource, _ptarget; |
|
144 Value _pstflow; |
|
132 |
145 |
133 // Result maps |
146 // Result maps |
134 FlowMap *_flow_map; |
147 FlowMap *_flow_map; |
135 PotentialMap *_potential_map; |
148 PotentialMap *_potential_map; |
136 bool _local_flow; |
149 bool _local_flow; |
137 bool _local_potential; |
150 bool _local_potential; |
138 |
151 |
139 // The number of nodes and arcs in the original graph |
152 // Data structures for storing the digraph |
140 int _node_num; |
|
141 int _arc_num; |
|
142 |
|
143 // Data structures for storing the graph |
|
144 IntNodeMap _node_id; |
153 IntNodeMap _node_id; |
145 ArcVector _arc_ref; |
154 ArcVector _arc_ref; |
146 IntVector _source; |
155 IntVector _source; |
147 IntVector _target; |
156 IntVector _target; |
148 |
157 |
149 // Node and arc maps |
158 // Node and arc data |
150 CapacityVector _cap; |
159 ValueVector _cap; |
151 CostVector _cost; |
160 ValueVector _cost; |
152 CostVector _supply; |
161 ValueVector _supply; |
153 CapacityVector _flow; |
162 ValueVector _flow; |
154 CostVector _pi; |
163 ValueVector _pi; |
155 |
164 |
156 // Data for storing the spanning tree structure |
165 // Data for storing the spanning tree structure |
157 IntVector _parent; |
166 IntVector _parent; |
158 IntVector _pred; |
167 IntVector _pred; |
159 IntVector _thread; |
168 IntVector _thread; |
167 |
176 |
168 // Temporary data used in the current pivot iteration |
177 // Temporary data used in the current pivot iteration |
169 int in_arc, join, u_in, v_in, u_out, v_out; |
178 int in_arc, join, u_in, v_in, u_out, v_out; |
170 int first, second, right, last; |
179 int first, second, right, last; |
171 int stem, par_stem, new_stem; |
180 int stem, par_stem, new_stem; |
172 Capacity delta; |
181 Value delta; |
173 |
182 |
174 private: |
183 private: |
175 |
184 |
176 /// \brief Implementation of the "First Eligible" pivot rule for the |
185 // Implementation of the First Eligible pivot rule |
177 /// \ref NetworkSimplex "network simplex" algorithm. |
|
178 /// |
|
179 /// This class implements the "First Eligible" pivot rule |
|
180 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
181 /// |
|
182 /// For more information see \ref NetworkSimplex::run(). |
|
183 class FirstEligiblePivotRule |
186 class FirstEligiblePivotRule |
184 { |
187 { |
185 private: |
188 private: |
186 |
189 |
187 // References to the NetworkSimplex class |
190 // References to the NetworkSimplex class |
188 const IntVector &_source; |
191 const IntVector &_source; |
189 const IntVector &_target; |
192 const IntVector &_target; |
190 const CostVector &_cost; |
193 const ValueVector &_cost; |
191 const IntVector &_state; |
194 const IntVector &_state; |
192 const CostVector &_pi; |
195 const ValueVector &_pi; |
193 int &_in_arc; |
196 int &_in_arc; |
194 int _arc_num; |
197 int _arc_num; |
195 |
198 |
196 // Pivot rule data |
199 // Pivot rule data |
197 int _next_arc; |
200 int _next_arc; |
198 |
201 |
199 public: |
202 public: |
200 |
203 |
201 /// Constructor |
204 // Constructor |
202 FirstEligiblePivotRule(NetworkSimplex &ns) : |
205 FirstEligiblePivotRule(NetworkSimplex &ns) : |
203 _source(ns._source), _target(ns._target), |
206 _source(ns._source), _target(ns._target), |
204 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
207 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
205 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
208 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
206 {} |
209 {} |
207 |
210 |
208 /// Find next entering arc |
211 // Find next entering arc |
209 bool findEnteringArc() { |
212 bool findEnteringArc() { |
210 Cost c; |
213 Value c; |
211 for (int e = _next_arc; e < _arc_num; ++e) { |
214 for (int e = _next_arc; e < _arc_num; ++e) { |
212 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
215 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
213 if (c < 0) { |
216 if (c < 0) { |
214 _in_arc = e; |
217 _in_arc = e; |
215 _next_arc = e + 1; |
218 _next_arc = e + 1; |
228 } |
231 } |
229 |
232 |
230 }; //class FirstEligiblePivotRule |
233 }; //class FirstEligiblePivotRule |
231 |
234 |
232 |
235 |
233 /// \brief Implementation of the "Best Eligible" pivot rule for the |
236 // Implementation of the Best Eligible pivot rule |
234 /// \ref NetworkSimplex "network simplex" algorithm. |
|
235 /// |
|
236 /// This class implements the "Best Eligible" pivot rule |
|
237 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
238 /// |
|
239 /// For more information see \ref NetworkSimplex::run(). |
|
240 class BestEligiblePivotRule |
237 class BestEligiblePivotRule |
241 { |
238 { |
242 private: |
239 private: |
243 |
240 |
244 // References to the NetworkSimplex class |
241 // References to the NetworkSimplex class |
245 const IntVector &_source; |
242 const IntVector &_source; |
246 const IntVector &_target; |
243 const IntVector &_target; |
247 const CostVector &_cost; |
244 const ValueVector &_cost; |
248 const IntVector &_state; |
245 const IntVector &_state; |
249 const CostVector &_pi; |
246 const ValueVector &_pi; |
250 int &_in_arc; |
247 int &_in_arc; |
251 int _arc_num; |
248 int _arc_num; |
252 |
249 |
253 public: |
250 public: |
254 |
251 |
255 /// Constructor |
252 // Constructor |
256 BestEligiblePivotRule(NetworkSimplex &ns) : |
253 BestEligiblePivotRule(NetworkSimplex &ns) : |
257 _source(ns._source), _target(ns._target), |
254 _source(ns._source), _target(ns._target), |
258 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
255 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
259 _in_arc(ns.in_arc), _arc_num(ns._arc_num) |
256 _in_arc(ns.in_arc), _arc_num(ns._arc_num) |
260 {} |
257 {} |
261 |
258 |
262 /// Find next entering arc |
259 // Find next entering arc |
263 bool findEnteringArc() { |
260 bool findEnteringArc() { |
264 Cost c, min = 0; |
261 Value c, min = 0; |
265 for (int e = 0; e < _arc_num; ++e) { |
262 for (int e = 0; e < _arc_num; ++e) { |
266 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
263 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
267 if (c < min) { |
264 if (c < min) { |
268 min = c; |
265 min = c; |
269 _in_arc = e; |
266 _in_arc = e; |
273 } |
270 } |
274 |
271 |
275 }; //class BestEligiblePivotRule |
272 }; //class BestEligiblePivotRule |
276 |
273 |
277 |
274 |
278 /// \brief Implementation of the "Block Search" pivot rule for the |
275 // Implementation of the Block Search pivot rule |
279 /// \ref NetworkSimplex "network simplex" algorithm. |
|
280 /// |
|
281 /// This class implements the "Block Search" pivot rule |
|
282 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
283 /// |
|
284 /// For more information see \ref NetworkSimplex::run(). |
|
285 class BlockSearchPivotRule |
276 class BlockSearchPivotRule |
286 { |
277 { |
287 private: |
278 private: |
288 |
279 |
289 // References to the NetworkSimplex class |
280 // References to the NetworkSimplex class |
290 const IntVector &_source; |
281 const IntVector &_source; |
291 const IntVector &_target; |
282 const IntVector &_target; |
292 const CostVector &_cost; |
283 const ValueVector &_cost; |
293 const IntVector &_state; |
284 const IntVector &_state; |
294 const CostVector &_pi; |
285 const ValueVector &_pi; |
295 int &_in_arc; |
286 int &_in_arc; |
296 int _arc_num; |
287 int _arc_num; |
297 |
288 |
298 // Pivot rule data |
289 // Pivot rule data |
299 int _block_size; |
290 int _block_size; |
300 int _next_arc; |
291 int _next_arc; |
301 |
292 |
302 public: |
293 public: |
303 |
294 |
304 /// Constructor |
295 // Constructor |
305 BlockSearchPivotRule(NetworkSimplex &ns) : |
296 BlockSearchPivotRule(NetworkSimplex &ns) : |
306 _source(ns._source), _target(ns._target), |
297 _source(ns._source), _target(ns._target), |
307 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
298 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
308 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
299 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
309 { |
300 { |
313 |
304 |
314 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)), |
305 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)), |
315 MIN_BLOCK_SIZE ); |
306 MIN_BLOCK_SIZE ); |
316 } |
307 } |
317 |
308 |
318 /// Find next entering arc |
309 // Find next entering arc |
319 bool findEnteringArc() { |
310 bool findEnteringArc() { |
320 Cost c, min = 0; |
311 Value c, min = 0; |
321 int cnt = _block_size; |
312 int cnt = _block_size; |
322 int e, min_arc = _next_arc; |
313 int e, min_arc = _next_arc; |
323 for (e = _next_arc; e < _arc_num; ++e) { |
314 for (e = _next_arc; e < _arc_num; ++e) { |
324 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
315 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
325 if (c < min) { |
316 if (c < min) { |
351 } |
342 } |
352 |
343 |
353 }; //class BlockSearchPivotRule |
344 }; //class BlockSearchPivotRule |
354 |
345 |
355 |
346 |
356 /// \brief Implementation of the "Candidate List" pivot rule for the |
347 // Implementation of the Candidate List pivot rule |
357 /// \ref NetworkSimplex "network simplex" algorithm. |
|
358 /// |
|
359 /// This class implements the "Candidate List" pivot rule |
|
360 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
361 /// |
|
362 /// For more information see \ref NetworkSimplex::run(). |
|
363 class CandidateListPivotRule |
348 class CandidateListPivotRule |
364 { |
349 { |
365 private: |
350 private: |
366 |
351 |
367 // References to the NetworkSimplex class |
352 // References to the NetworkSimplex class |
368 const IntVector &_source; |
353 const IntVector &_source; |
369 const IntVector &_target; |
354 const IntVector &_target; |
370 const CostVector &_cost; |
355 const ValueVector &_cost; |
371 const IntVector &_state; |
356 const IntVector &_state; |
372 const CostVector &_pi; |
357 const ValueVector &_pi; |
373 int &_in_arc; |
358 int &_in_arc; |
374 int _arc_num; |
359 int _arc_num; |
375 |
360 |
376 // Pivot rule data |
361 // Pivot rule data |
377 IntVector _candidates; |
362 IntVector _candidates; |
401 _candidates.resize(_list_length); |
386 _candidates.resize(_list_length); |
402 } |
387 } |
403 |
388 |
404 /// Find next entering arc |
389 /// Find next entering arc |
405 bool findEnteringArc() { |
390 bool findEnteringArc() { |
406 Cost min, c; |
391 Value min, c; |
407 int e, min_arc = _next_arc; |
392 int e, min_arc = _next_arc; |
408 if (_curr_length > 0 && _minor_count < _minor_limit) { |
393 if (_curr_length > 0 && _minor_count < _minor_limit) { |
409 // Minor iteration: select the best eligible arc from the |
394 // Minor iteration: select the best eligible arc from the |
410 // current candidate list |
395 // current candidate list |
411 ++_minor_count; |
396 ++_minor_count; |
462 } |
447 } |
463 |
448 |
464 }; //class CandidateListPivotRule |
449 }; //class CandidateListPivotRule |
465 |
450 |
466 |
451 |
467 /// \brief Implementation of the "Altering Candidate List" pivot rule |
452 // Implementation of the Altering Candidate List pivot rule |
468 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
469 /// |
|
470 /// This class implements the "Altering Candidate List" pivot rule |
|
471 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
472 /// |
|
473 /// For more information see \ref NetworkSimplex::run(). |
|
474 class AlteringListPivotRule |
453 class AlteringListPivotRule |
475 { |
454 { |
476 private: |
455 private: |
477 |
456 |
478 // References to the NetworkSimplex class |
457 // References to the NetworkSimplex class |
479 const IntVector &_source; |
458 const IntVector &_source; |
480 const IntVector &_target; |
459 const IntVector &_target; |
481 const CostVector &_cost; |
460 const ValueVector &_cost; |
482 const IntVector &_state; |
461 const IntVector &_state; |
483 const CostVector &_pi; |
462 const ValueVector &_pi; |
484 int &_in_arc; |
463 int &_in_arc; |
485 int _arc_num; |
464 int _arc_num; |
486 |
465 |
487 // Pivot rule data |
466 // Pivot rule data |
488 int _block_size, _head_length, _curr_length; |
467 int _block_size, _head_length, _curr_length; |
489 int _next_arc; |
468 int _next_arc; |
490 IntVector _candidates; |
469 IntVector _candidates; |
491 CostVector _cand_cost; |
470 ValueVector _cand_cost; |
492 |
471 |
493 // Functor class to compare arcs during sort of the candidate list |
472 // Functor class to compare arcs during sort of the candidate list |
494 class SortFunc |
473 class SortFunc |
495 { |
474 { |
496 private: |
475 private: |
497 const CostVector &_map; |
476 const ValueVector &_map; |
498 public: |
477 public: |
499 SortFunc(const CostVector &map) : _map(map) {} |
478 SortFunc(const ValueVector &map) : _map(map) {} |
500 bool operator()(int left, int right) { |
479 bool operator()(int left, int right) { |
501 return _map[left] > _map[right]; |
480 return _map[left] > _map[right]; |
502 } |
481 } |
503 }; |
482 }; |
504 |
483 |
505 SortFunc _sort_func; |
484 SortFunc _sort_func; |
506 |
485 |
507 public: |
486 public: |
508 |
487 |
509 /// Constructor |
488 // Constructor |
510 AlteringListPivotRule(NetworkSimplex &ns) : |
489 AlteringListPivotRule(NetworkSimplex &ns) : |
511 _source(ns._source), _target(ns._target), |
490 _source(ns._source), _target(ns._target), |
512 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
491 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
513 _in_arc(ns.in_arc), _arc_num(ns._arc_num), |
492 _in_arc(ns.in_arc), _arc_num(ns._arc_num), |
514 _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost) |
493 _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost) |
525 MIN_HEAD_LENGTH ); |
504 MIN_HEAD_LENGTH ); |
526 _candidates.resize(_head_length + _block_size); |
505 _candidates.resize(_head_length + _block_size); |
527 _curr_length = 0; |
506 _curr_length = 0; |
528 } |
507 } |
529 |
508 |
530 /// Find next entering arc |
509 // Find next entering arc |
531 bool findEnteringArc() { |
510 bool findEnteringArc() { |
532 // Check the current candidate list |
511 // Check the current candidate list |
533 int e; |
512 int e; |
534 for (int i = 0; i < _curr_length; ++i) { |
513 for (int i = 0; i < _curr_length; ++i) { |
535 e = _candidates[i]; |
514 e = _candidates[i]; |
590 |
569 |
591 }; //class AlteringListPivotRule |
570 }; //class AlteringListPivotRule |
592 |
571 |
593 public: |
572 public: |
594 |
573 |
595 /// \brief General constructor (with lower bounds). |
574 /// \brief Constructor. |
596 /// |
575 /// |
597 /// General constructor (with lower bounds). |
576 /// Constructor. |
598 /// |
577 /// |
599 /// \param graph The digraph the algorithm runs on. |
578 /// \param graph The digraph the algorithm runs on. |
600 /// \param lower The lower bounds of the arcs. |
579 NetworkSimplex(const GR& graph) : |
601 /// \param capacity The capacities (upper bounds) of the arcs. |
580 _graph(graph), |
602 /// \param cost The cost (length) values of the arcs. |
581 _plower(NULL), _pupper(NULL), _pcost(NULL), |
603 /// \param supply The supply values of the nodes (signed). |
582 _psupply(NULL), _pstsup(false), |
604 NetworkSimplex( const Digraph &graph, |
|
605 const LowerMap &lower, |
|
606 const CapacityMap &capacity, |
|
607 const CostMap &cost, |
|
608 const SupplyMap &supply ) : |
|
609 _graph(graph), _orig_lower(&lower), _orig_cap(capacity), |
|
610 _orig_cost(cost), _orig_supply(&supply), |
|
611 _flow_map(NULL), _potential_map(NULL), |
583 _flow_map(NULL), _potential_map(NULL), |
612 _local_flow(false), _local_potential(false), |
584 _local_flow(false), _local_potential(false), |
613 _node_id(graph) |
585 _node_id(graph) |
614 {} |
586 { |
615 |
587 LEMON_ASSERT(std::numeric_limits<Value>::is_integer && |
616 /// \brief General constructor (without lower bounds). |
588 std::numeric_limits<Value>::is_signed, |
617 /// |
589 "The value type of NetworkSimplex must be a signed integer"); |
618 /// General constructor (without lower bounds). |
590 } |
619 /// |
|
620 /// \param graph The digraph the algorithm runs on. |
|
621 /// \param capacity The capacities (upper bounds) of the arcs. |
|
622 /// \param cost The cost (length) values of the arcs. |
|
623 /// \param supply The supply values of the nodes (signed). |
|
624 NetworkSimplex( const Digraph &graph, |
|
625 const CapacityMap &capacity, |
|
626 const CostMap &cost, |
|
627 const SupplyMap &supply ) : |
|
628 _graph(graph), _orig_lower(NULL), _orig_cap(capacity), |
|
629 _orig_cost(cost), _orig_supply(&supply), |
|
630 _flow_map(NULL), _potential_map(NULL), |
|
631 _local_flow(false), _local_potential(false), |
|
632 _node_id(graph) |
|
633 {} |
|
634 |
|
635 /// \brief Simple constructor (with lower bounds). |
|
636 /// |
|
637 /// Simple constructor (with lower bounds). |
|
638 /// |
|
639 /// \param graph The digraph the algorithm runs on. |
|
640 /// \param lower The lower bounds of the arcs. |
|
641 /// \param capacity The capacities (upper bounds) of the arcs. |
|
642 /// \param cost The cost (length) values of the arcs. |
|
643 /// \param s The source node. |
|
644 /// \param t The target node. |
|
645 /// \param flow_value The required amount of flow from node \c s |
|
646 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
|
647 NetworkSimplex( const Digraph &graph, |
|
648 const LowerMap &lower, |
|
649 const CapacityMap &capacity, |
|
650 const CostMap &cost, |
|
651 Node s, Node t, |
|
652 Capacity flow_value ) : |
|
653 _graph(graph), _orig_lower(&lower), _orig_cap(capacity), |
|
654 _orig_cost(cost), _orig_supply(NULL), |
|
655 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), |
|
656 _flow_map(NULL), _potential_map(NULL), |
|
657 _local_flow(false), _local_potential(false), |
|
658 _node_id(graph) |
|
659 {} |
|
660 |
|
661 /// \brief Simple constructor (without lower bounds). |
|
662 /// |
|
663 /// Simple constructor (without lower bounds). |
|
664 /// |
|
665 /// \param graph The digraph the algorithm runs on. |
|
666 /// \param capacity The capacities (upper bounds) of the arcs. |
|
667 /// \param cost The cost (length) values of the arcs. |
|
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 Digraph &graph, |
|
673 const CapacityMap &capacity, |
|
674 const CostMap &cost, |
|
675 Node s, Node t, |
|
676 Capacity flow_value ) : |
|
677 _graph(graph), _orig_lower(NULL), _orig_cap(capacity), |
|
678 _orig_cost(cost), _orig_supply(NULL), |
|
679 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), |
|
680 _flow_map(NULL), _potential_map(NULL), |
|
681 _local_flow(false), _local_potential(false), |
|
682 _node_id(graph) |
|
683 {} |
|
684 |
591 |
685 /// Destructor. |
592 /// Destructor. |
686 ~NetworkSimplex() { |
593 ~NetworkSimplex() { |
687 if (_local_flow) delete _flow_map; |
594 if (_local_flow) delete _flow_map; |
688 if (_local_potential) delete _potential_map; |
595 if (_local_potential) delete _potential_map; |
689 } |
596 } |
690 |
597 |
598 /// \brief Set the lower bounds on the arcs. |
|
599 /// |
|
600 /// This function sets the lower bounds on the arcs. |
|
601 /// If neither this function nor \ref boundMaps() is used before |
|
602 /// calling \ref run(), the lower bounds will be set to zero |
|
603 /// on all arcs. |
|
604 /// |
|
605 /// \param map An arc map storing the lower bounds. |
|
606 /// Its \c Value type must be convertible to the \c Value type |
|
607 /// of the algorithm. |
|
608 /// |
|
609 /// \return <tt>(*this)</tt> |
|
610 template <typename LOWER> |
|
611 NetworkSimplex& lowerMap(const LOWER& map) { |
|
612 delete _plower; |
|
613 _plower = new ValueArcMap(_graph); |
|
614 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
615 (*_plower)[a] = map[a]; |
|
616 } |
|
617 return *this; |
|
618 } |
|
619 |
|
620 /// \brief Set the upper bounds (capacities) on the arcs. |
|
621 /// |
|
622 /// This function sets the upper bounds (capacities) on the arcs. |
|
623 /// If none of the functions \ref upperMap(), \ref capacityMap() |
|
624 /// and \ref boundMaps() is used before calling \ref run(), |
|
625 /// the upper bounds (capacities) will be set to |
|
626 /// \c std::numeric_limits<Value>::max() on all arcs. |
|
627 /// |
|
628 /// \param map An arc map storing the upper bounds. |
|
629 /// Its \c Value type must be convertible to the \c Value type |
|
630 /// of the algorithm. |
|
631 /// |
|
632 /// \return <tt>(*this)</tt> |
|
633 template<typename UPPER> |
|
634 NetworkSimplex& upperMap(const UPPER& map) { |
|
635 delete _pupper; |
|
636 _pupper = new ValueArcMap(_graph); |
|
637 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
638 (*_pupper)[a] = map[a]; |
|
639 } |
|
640 return *this; |
|
641 } |
|
642 |
|
643 /// \brief Set the upper bounds (capacities) on the arcs. |
|
644 /// |
|
645 /// This function sets the upper bounds (capacities) on the arcs. |
|
646 /// It is just an alias for \ref upperMap(). |
|
647 /// |
|
648 /// \return <tt>(*this)</tt> |
|
649 template<typename CAP> |
|
650 NetworkSimplex& capacityMap(const CAP& map) { |
|
651 return upperMap(map); |
|
652 } |
|
653 |
|
654 /// \brief Set the lower and upper bounds on the arcs. |
|
655 /// |
|
656 /// This function sets the lower and upper bounds on the arcs. |
|
657 /// If neither this function nor \ref lowerMap() is used before |
|
658 /// calling \ref run(), the lower bounds will be set to zero |
|
659 /// on all arcs. |
|
660 /// If none of the functions \ref upperMap(), \ref capacityMap() |
|
661 /// and \ref boundMaps() is used before calling \ref run(), |
|
662 /// the upper bounds (capacities) will be set to |
|
663 /// \c std::numeric_limits<Value>::max() on all arcs. |
|
664 /// |
|
665 /// \param lower An arc map storing the lower bounds. |
|
666 /// \param upper An arc map storing the upper bounds. |
|
667 /// |
|
668 /// The \c Value type of the maps must be convertible to the |
|
669 /// \c Value type of the algorithm. |
|
670 /// |
|
671 /// \note This function is just a shortcut of calling \ref lowerMap() |
|
672 /// and \ref upperMap() separately. |
|
673 /// |
|
674 /// \return <tt>(*this)</tt> |
|
675 template <typename LOWER, typename UPPER> |
|
676 NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) { |
|
677 return lowerMap(lower).upperMap(upper); |
|
678 } |
|
679 |
|
680 /// \brief Set the costs of the arcs. |
|
681 /// |
|
682 /// This function sets the costs of the arcs. |
|
683 /// If it is not used before calling \ref run(), the costs |
|
684 /// will be set to \c 1 on all arcs. |
|
685 /// |
|
686 /// \param map An arc map storing the costs. |
|
687 /// Its \c Value type must be convertible to the \c Value type |
|
688 /// of the algorithm. |
|
689 /// |
|
690 /// \return <tt>(*this)</tt> |
|
691 template<typename COST> |
|
692 NetworkSimplex& costMap(const COST& map) { |
|
693 delete _pcost; |
|
694 _pcost = new ValueArcMap(_graph); |
|
695 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
696 (*_pcost)[a] = map[a]; |
|
697 } |
|
698 return *this; |
|
699 } |
|
700 |
|
701 /// \brief Set the supply values of the nodes. |
|
702 /// |
|
703 /// This function sets the supply values of the nodes. |
|
704 /// If neither this function nor \ref stSupply() is used before |
|
705 /// calling \ref run(), the supply of each node will be set to zero. |
|
706 /// (It makes sense only if non-zero lower bounds are given.) |
|
707 /// |
|
708 /// \param map A node map storing the supply values. |
|
709 /// Its \c Value type must be convertible to the \c Value type |
|
710 /// of the algorithm. |
|
711 /// |
|
712 /// \return <tt>(*this)</tt> |
|
713 template<typename SUP> |
|
714 NetworkSimplex& supplyMap(const SUP& map) { |
|
715 delete _psupply; |
|
716 _pstsup = false; |
|
717 _psupply = new ValueNodeMap(_graph); |
|
718 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
719 (*_psupply)[n] = map[n]; |
|
720 } |
|
721 return *this; |
|
722 } |
|
723 |
|
724 /// \brief Set single source and target nodes and a supply value. |
|
725 /// |
|
726 /// This function sets a single source node and a single target node |
|
727 /// and the required flow value. |
|
728 /// If neither this function nor \ref supplyMap() is used before |
|
729 /// calling \ref run(), the supply of each node will be set to zero. |
|
730 /// (It makes sense only if non-zero lower bounds are given.) |
|
731 /// |
|
732 /// \param s The source node. |
|
733 /// \param t The target node. |
|
734 /// \param k The required amount of flow from node \c s to node \c t |
|
735 /// (i.e. the supply of \c s and the demand of \c t). |
|
736 /// |
|
737 /// \return <tt>(*this)</tt> |
|
738 NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { |
|
739 delete _psupply; |
|
740 _psupply = NULL; |
|
741 _pstsup = true; |
|
742 _psource = s; |
|
743 _ptarget = t; |
|
744 _pstflow = k; |
|
745 return *this; |
|
746 } |
|
747 |
|
691 /// \brief Set the flow map. |
748 /// \brief Set the flow map. |
692 /// |
749 /// |
693 /// This function sets the flow map. |
750 /// This function sets the flow map. |
751 /// If it is not used before calling \ref run(), an instance will |
|
752 /// be allocated automatically. The destructor deallocates this |
|
753 /// automatically allocated map, of course. |
|
694 /// |
754 /// |
695 /// \return <tt>(*this)</tt> |
755 /// \return <tt>(*this)</tt> |
696 NetworkSimplex& flowMap(FlowMap &map) { |
756 NetworkSimplex& flowMap(FlowMap& map) { |
697 if (_local_flow) { |
757 if (_local_flow) { |
698 delete _flow_map; |
758 delete _flow_map; |
699 _local_flow = false; |
759 _local_flow = false; |
700 } |
760 } |
701 _flow_map = ↦ |
761 _flow_map = ↦ |
702 return *this; |
762 return *this; |
703 } |
763 } |
704 |
764 |
705 /// \brief Set the potential map. |
765 /// \brief Set the potential map. |
706 /// |
766 /// |
707 /// This function sets the potential map. |
767 /// This function sets the potential map, which is used for storing |
768 /// the dual solution. |
|
769 /// If it is not used before calling \ref run(), an instance will |
|
770 /// be allocated automatically. The destructor deallocates this |
|
771 /// automatically allocated map, of course. |
|
708 /// |
772 /// |
709 /// \return <tt>(*this)</tt> |
773 /// \return <tt>(*this)</tt> |
710 NetworkSimplex& potentialMap(PotentialMap &map) { |
774 NetworkSimplex& potentialMap(PotentialMap& map) { |
711 if (_local_potential) { |
775 if (_local_potential) { |
712 delete _potential_map; |
776 delete _potential_map; |
713 _local_potential = false; |
777 _local_potential = false; |
714 } |
778 } |
715 _potential_map = ↦ |
779 _potential_map = ↦ |
716 return *this; |
780 return *this; |
717 } |
781 } |
718 |
782 |
719 /// \name Execution control |
783 /// \name Execution Control |
720 /// The algorithm can be executed using the |
784 /// The algorithm can be executed using \ref run(). |
721 /// \ref lemon::NetworkSimplex::run() "run()" function. |
785 |
722 /// @{ |
786 /// @{ |
723 |
787 |
724 /// \brief Run the algorithm. |
788 /// \brief Run the algorithm. |
725 /// |
789 /// |
726 /// This function runs the algorithm. |
790 /// This function runs the algorithm. |
727 /// |
791 /// The paramters can be specified using \ref lowerMap(), |
728 /// \param pivot_rule The pivot rule that is used during the |
792 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), |
729 /// algorithm. |
793 /// \ref costMap(), \ref supplyMap() and \ref stSupply() |
730 /// |
794 /// functions. For example, |
731 /// The available pivot rules: |
795 /// \code |
732 /// |
796 /// NetworkSimplex<ListDigraph> ns(graph); |
733 /// - FIRST_ELIGIBLE_PIVOT The next eligible arc is selected in |
797 /// ns.boundMaps(lower, upper).costMap(cost) |
734 /// a wraparound fashion in every iteration |
798 /// .supplyMap(sup).run(); |
735 /// (\ref FirstEligiblePivotRule). |
799 /// \endcode |
736 /// |
800 /// |
737 /// - BEST_ELIGIBLE_PIVOT The best eligible arc is selected in |
801 /// \param pivot_rule The pivot rule that will be used during the |
738 /// every iteration (\ref BestEligiblePivotRule). |
802 /// algorithm. For more information see \ref PivotRule. |
739 /// |
|
740 /// - BLOCK_SEARCH_PIVOT A specified number of arcs are examined in |
|
741 /// every iteration in a wraparound fashion and the best eligible |
|
742 /// arc is selected from this block (\ref BlockSearchPivotRule). |
|
743 /// |
|
744 /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is |
|
745 /// built from eligible arcs in a wraparound fashion and in the |
|
746 /// following minor iterations the best eligible arc is selected |
|
747 /// from this list (\ref CandidateListPivotRule). |
|
748 /// |
|
749 /// - ALTERING_LIST_PIVOT It is a modified version of the |
|
750 /// "Candidate List" pivot rule. It keeps only the several best |
|
751 /// eligible arcs from the former candidate list and extends this |
|
752 /// list in every iteration (\ref AlteringListPivotRule). |
|
753 /// |
|
754 /// According to our comprehensive benchmark tests the "Block Search" |
|
755 /// pivot rule proved to be the fastest and the most robust on |
|
756 /// various test inputs. Thus it is the default option. |
|
757 /// |
803 /// |
758 /// \return \c true if a feasible flow can be found. |
804 /// \return \c true if a feasible flow can be found. |
759 bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) { |
805 bool run(PivotRule pivot_rule = BLOCK_SEARCH) { |
760 return init() && start(pivot_rule); |
806 return init() && start(pivot_rule); |
761 } |
807 } |
762 |
808 |
763 /// @} |
809 /// @} |
764 |
810 |
765 /// \name Query Functions |
811 /// \name Query Functions |
766 /// The results of the algorithm can be obtained using these |
812 /// The results of the algorithm can be obtained using these |
767 /// functions.\n |
813 /// functions.\n |
768 /// \ref lemon::NetworkSimplex::run() "run()" must be called before |
814 /// The \ref run() function must be called before using them. |
769 /// using them. |
815 |
770 /// @{ |
816 /// @{ |
817 |
|
818 /// \brief Return the total cost of the found flow. |
|
819 /// |
|
820 /// This function returns the total cost of the found flow. |
|
821 /// The complexity of the function is \f$ O(e) \f$. |
|
822 /// |
|
823 /// \note The return type of the function can be specified as a |
|
824 /// template parameter. For example, |
|
825 /// \code |
|
826 /// ns.totalCost<double>(); |
|
827 /// \endcode |
|
828 /// It is useful if the total cost cannot be stored in the \c Value |
|
829 /// type of the algorithm, which is the default return type of the |
|
830 /// function. |
|
831 /// |
|
832 /// \pre \ref run() must be called before using this function. |
|
833 template <typename Num> |
|
834 Num totalCost() const { |
|
835 Num c = 0; |
|
836 if (_pcost) { |
|
837 for (ArcIt e(_graph); e != INVALID; ++e) |
|
838 c += (*_flow_map)[e] * (*_pcost)[e]; |
|
839 } else { |
|
840 for (ArcIt e(_graph); e != INVALID; ++e) |
|
841 c += (*_flow_map)[e]; |
|
842 } |
|
843 return c; |
|
844 } |
|
845 |
|
846 #ifndef DOXYGEN |
|
847 Value totalCost() const { |
|
848 return totalCost<Value>(); |
|
849 } |
|
850 #endif |
|
851 |
|
852 /// \brief Return the flow on the given arc. |
|
853 /// |
|
854 /// This function returns the flow on the given arc. |
|
855 /// |
|
856 /// \pre \ref run() must be called before using this function. |
|
857 Value flow(const Arc& a) const { |
|
858 return (*_flow_map)[a]; |
|
859 } |
|
771 |
860 |
772 /// \brief Return a const reference to the flow map. |
861 /// \brief Return a const reference to the flow map. |
773 /// |
862 /// |
774 /// This function returns a const reference to an arc map storing |
863 /// This function returns a const reference to an arc map storing |
775 /// the found flow. |
864 /// the found flow. |
777 /// \pre \ref run() must be called before using this function. |
866 /// \pre \ref run() must be called before using this function. |
778 const FlowMap& flowMap() const { |
867 const FlowMap& flowMap() const { |
779 return *_flow_map; |
868 return *_flow_map; |
780 } |
869 } |
781 |
870 |
871 /// \brief Return the potential (dual value) of the given node. |
|
872 /// |
|
873 /// This function returns the potential (dual value) of the |
|
874 /// given node. |
|
875 /// |
|
876 /// \pre \ref run() must be called before using this function. |
|
877 Value potential(const Node& n) const { |
|
878 return (*_potential_map)[n]; |
|
879 } |
|
880 |
|
782 /// \brief Return a const reference to the potential map |
881 /// \brief Return a const reference to the potential map |
783 /// (the dual solution). |
882 /// (the dual solution). |
784 /// |
883 /// |
785 /// This function returns a const reference to a node map storing |
884 /// This function returns a const reference to a node map storing |
786 /// the found potentials (the dual solution). |
885 /// the found potentials, which form the dual solution of the |
886 /// \ref min_cost_flow "minimum cost flow" problem. |
|
787 /// |
887 /// |
788 /// \pre \ref run() must be called before using this function. |
888 /// \pre \ref run() must be called before using this function. |
789 const PotentialMap& potentialMap() const { |
889 const PotentialMap& potentialMap() const { |
790 return *_potential_map; |
890 return *_potential_map; |
791 } |
|
792 |
|
793 /// \brief Return the flow on the given arc. |
|
794 /// |
|
795 /// This function returns the flow on the given arc. |
|
796 /// |
|
797 /// \pre \ref run() must be called before using this function. |
|
798 Capacity flow(const Arc& arc) const { |
|
799 return (*_flow_map)[arc]; |
|
800 } |
|
801 |
|
802 /// \brief Return the potential of the given node. |
|
803 /// |
|
804 /// This function returns the potential of the given node. |
|
805 /// |
|
806 /// \pre \ref run() must be called before using this function. |
|
807 Cost potential(const Node& node) const { |
|
808 return (*_potential_map)[node]; |
|
809 } |
|
810 |
|
811 /// \brief Return the total cost of the found flow. |
|
812 /// |
|
813 /// This function returns the total cost of the found flow. |
|
814 /// The complexity of the function is \f$ O(e) \f$. |
|
815 /// |
|
816 /// \pre \ref run() must be called before using this function. |
|
817 Cost totalCost() const { |
|
818 Cost c = 0; |
|
819 for (ArcIt e(_graph); e != INVALID; ++e) |
|
820 c += (*_flow_map)[e] * _orig_cost[e]; |
|
821 return c; |
|
822 } |
891 } |
823 |
892 |
824 /// @} |
893 /// @} |
825 |
894 |
826 private: |
895 private: |
840 // Initialize vectors |
909 // Initialize vectors |
841 _node_num = countNodes(_graph); |
910 _node_num = countNodes(_graph); |
842 _arc_num = countArcs(_graph); |
911 _arc_num = countArcs(_graph); |
843 int all_node_num = _node_num + 1; |
912 int all_node_num = _node_num + 1; |
844 int all_arc_num = _arc_num + _node_num; |
913 int all_arc_num = _arc_num + _node_num; |
914 if (_node_num == 0) return false; |
|
845 |
915 |
846 _arc_ref.resize(_arc_num); |
916 _arc_ref.resize(_arc_num); |
847 _source.resize(all_arc_num); |
917 _source.resize(all_arc_num); |
848 _target.resize(all_arc_num); |
918 _target.resize(all_arc_num); |
849 |
919 |
862 _last_succ.resize(all_node_num); |
932 _last_succ.resize(all_node_num); |
863 _state.resize(all_arc_num, STATE_LOWER); |
933 _state.resize(all_arc_num, STATE_LOWER); |
864 |
934 |
865 // Initialize node related data |
935 // Initialize node related data |
866 bool valid_supply = true; |
936 bool valid_supply = true; |
867 if (_orig_supply) { |
937 if (!_pstsup && !_psupply) { |
868 Supply sum = 0; |
938 _pstsup = true; |
939 _psource = _ptarget = NodeIt(_graph); |
|
940 _pstflow = 0; |
|
941 } |
|
942 if (_psupply) { |
|
943 Value sum = 0; |
|
869 int i = 0; |
944 int i = 0; |
870 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
945 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
871 _node_id[n] = i; |
946 _node_id[n] = i; |
872 _supply[i] = (*_orig_supply)[n]; |
947 _supply[i] = (*_psupply)[n]; |
873 sum += _supply[i]; |
948 sum += _supply[i]; |
874 } |
949 } |
875 valid_supply = (sum == 0); |
950 valid_supply = (sum == 0); |
876 } else { |
951 } else { |
877 int i = 0; |
952 int i = 0; |
878 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
953 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
879 _node_id[n] = i; |
954 _node_id[n] = i; |
880 _supply[i] = 0; |
955 _supply[i] = 0; |
881 } |
956 } |
882 _supply[_node_id[_orig_source]] = _orig_flow_value; |
957 _supply[_node_id[_psource]] = _pstflow; |
883 _supply[_node_id[_orig_target]] = -_orig_flow_value; |
958 _supply[_node_id[_ptarget]] = -_pstflow; |
884 } |
959 } |
885 if (!valid_supply) return false; |
960 if (!valid_supply) return false; |
886 |
961 |
887 // Set data for the artificial root node |
962 // Set data for the artificial root node |
888 _root = _node_num; |
963 _root = _node_num; |
902 _arc_ref[i] = e; |
977 _arc_ref[i] = e; |
903 if ((i += k) >= _arc_num) i = (i % k) + 1; |
978 if ((i += k) >= _arc_num) i = (i % k) + 1; |
904 } |
979 } |
905 |
980 |
906 // Initialize arc maps |
981 // Initialize arc maps |
907 for (int i = 0; i != _arc_num; ++i) { |
982 if (_pupper && _pcost) { |
908 Arc e = _arc_ref[i]; |
983 for (int i = 0; i != _arc_num; ++i) { |
909 _source[i] = _node_id[_graph.source(e)]; |
984 Arc e = _arc_ref[i]; |
910 _target[i] = _node_id[_graph.target(e)]; |
985 _source[i] = _node_id[_graph.source(e)]; |
911 _cost[i] = _orig_cost[e]; |
986 _target[i] = _node_id[_graph.target(e)]; |
912 _cap[i] = _orig_cap[e]; |
987 _cap[i] = (*_pupper)[e]; |
988 _cost[i] = (*_pcost)[e]; |
|
989 } |
|
990 } else { |
|
991 for (int i = 0; i != _arc_num; ++i) { |
|
992 Arc e = _arc_ref[i]; |
|
993 _source[i] = _node_id[_graph.source(e)]; |
|
994 _target[i] = _node_id[_graph.target(e)]; |
|
995 } |
|
996 if (_pupper) { |
|
997 for (int i = 0; i != _arc_num; ++i) |
|
998 _cap[i] = (*_pupper)[_arc_ref[i]]; |
|
999 } else { |
|
1000 Value val = std::numeric_limits<Value>::max(); |
|
1001 for (int i = 0; i != _arc_num; ++i) |
|
1002 _cap[i] = val; |
|
1003 } |
|
1004 if (_pcost) { |
|
1005 for (int i = 0; i != _arc_num; ++i) |
|
1006 _cost[i] = (*_pcost)[_arc_ref[i]]; |
|
1007 } else { |
|
1008 for (int i = 0; i != _arc_num; ++i) |
|
1009 _cost[i] = 1; |
|
1010 } |
|
913 } |
1011 } |
914 |
1012 |
915 // Remove non-zero lower bounds |
1013 // Remove non-zero lower bounds |
916 if (_orig_lower) { |
1014 if (_plower) { |
917 for (int i = 0; i != _arc_num; ++i) { |
1015 for (int i = 0; i != _arc_num; ++i) { |
918 Capacity c = (*_orig_lower)[_arc_ref[i]]; |
1016 Value c = (*_plower)[_arc_ref[i]]; |
919 if (c != 0) { |
1017 if (c != 0) { |
920 _cap[i] -= c; |
1018 _cap[i] -= c; |
921 _supply[_source[i]] -= c; |
1019 _supply[_source[i]] -= c; |
922 _supply[_target[i]] += c; |
1020 _supply[_target[i]] += c; |
923 } |
1021 } |
924 } |
1022 } |
925 } |
1023 } |
926 |
1024 |
927 // Add artificial arcs and initialize the spanning tree data structure |
1025 // Add artificial arcs and initialize the spanning tree data structure |
928 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
1026 Value max_cap = std::numeric_limits<Value>::max(); |
929 Capacity max_cap = std::numeric_limits<Capacity>::max(); |
1027 Value max_cost = std::numeric_limits<Value>::max() / 4; |
930 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
1028 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
931 _thread[u] = u + 1; |
1029 _thread[u] = u + 1; |
932 _rev_thread[u + 1] = u; |
1030 _rev_thread[u + 1] = u; |
933 _succ_num[u] = 1; |
1031 _succ_num[u] = 1; |
934 _last_succ[u] = u; |
1032 _last_succ[u] = u; |
977 first = _target[in_arc]; |
1075 first = _target[in_arc]; |
978 second = _source[in_arc]; |
1076 second = _source[in_arc]; |
979 } |
1077 } |
980 delta = _cap[in_arc]; |
1078 delta = _cap[in_arc]; |
981 int result = 0; |
1079 int result = 0; |
982 Capacity d; |
1080 Value d; |
983 int e; |
1081 int e; |
984 |
1082 |
985 // Search the cycle along the path form the first node to the root |
1083 // Search the cycle along the path form the first node to the root |
986 for (int u = first; u != join; u = _parent[u]) { |
1084 for (int u = first; u != join; u = _parent[u]) { |
987 e = _pred[u]; |
1085 e = _pred[u]; |
1015 |
1113 |
1016 // Change _flow and _state vectors |
1114 // Change _flow and _state vectors |
1017 void changeFlow(bool change) { |
1115 void changeFlow(bool change) { |
1018 // Augment along the cycle |
1116 // Augment along the cycle |
1019 if (delta > 0) { |
1117 if (delta > 0) { |
1020 Capacity val = _state[in_arc] * delta; |
1118 Value val = _state[in_arc] * delta; |
1021 _flow[in_arc] += val; |
1119 _flow[in_arc] += val; |
1022 for (int u = _source[in_arc]; u != join; u = _parent[u]) { |
1120 for (int u = _source[in_arc]; u != join; u = _parent[u]) { |
1023 _flow[_pred[u]] += _forward[u] ? -val : val; |
1121 _flow[_pred[u]] += _forward[u] ? -val : val; |
1024 } |
1122 } |
1025 for (int u = _target[in_arc]; u != join; u = _parent[u]) { |
1123 for (int u = _target[in_arc]; u != join; u = _parent[u]) { |
1156 } |
1254 } |
1157 } |
1255 } |
1158 |
1256 |
1159 // Update potentials |
1257 // Update potentials |
1160 void updatePotential() { |
1258 void updatePotential() { |
1161 Cost sigma = _forward[u_in] ? |
1259 Value sigma = _forward[u_in] ? |
1162 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1260 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1163 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1261 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1164 if (_succ_num[u_in] > _node_num / 2) { |
1262 if (_succ_num[u_in] > _node_num / 2) { |
1165 // Update in the upper subtree (which contains the root) |
1263 // Update in the upper subtree (which contains the root) |
1166 int before = _rev_thread[u_in]; |
1264 int before = _rev_thread[u_in]; |
1179 } |
1277 } |
1180 } |
1278 } |
1181 } |
1279 } |
1182 |
1280 |
1183 // Execute the algorithm |
1281 // Execute the algorithm |
1184 bool start(PivotRuleEnum pivot_rule) { |
1282 bool start(PivotRule pivot_rule) { |
1185 // Select the pivot rule implementation |
1283 // Select the pivot rule implementation |
1186 switch (pivot_rule) { |
1284 switch (pivot_rule) { |
1187 case FIRST_ELIGIBLE_PIVOT: |
1285 case FIRST_ELIGIBLE: |
1188 return start<FirstEligiblePivotRule>(); |
1286 return start<FirstEligiblePivotRule>(); |
1189 case BEST_ELIGIBLE_PIVOT: |
1287 case BEST_ELIGIBLE: |
1190 return start<BestEligiblePivotRule>(); |
1288 return start<BestEligiblePivotRule>(); |
1191 case BLOCK_SEARCH_PIVOT: |
1289 case BLOCK_SEARCH: |
1192 return start<BlockSearchPivotRule>(); |
1290 return start<BlockSearchPivotRule>(); |
1193 case CANDIDATE_LIST_PIVOT: |
1291 case CANDIDATE_LIST: |
1194 return start<CandidateListPivotRule>(); |
1292 return start<CandidateListPivotRule>(); |
1195 case ALTERING_LIST_PIVOT: |
1293 case ALTERING_LIST: |
1196 return start<AlteringListPivotRule>(); |
1294 return start<AlteringListPivotRule>(); |
1197 } |
1295 } |
1198 return false; |
1296 return false; |
1199 } |
1297 } |
1200 |
1298 |
1201 template<class PivotRuleImplementation> |
1299 template <typename PivotRuleImpl> |
1202 bool start() { |
1300 bool start() { |
1203 PivotRuleImplementation pivot(*this); |
1301 PivotRuleImpl pivot(*this); |
1204 |
1302 |
1205 // Execute the network simplex algorithm |
1303 // Execute the Network Simplex algorithm |
1206 while (pivot.findEnteringArc()) { |
1304 while (pivot.findEnteringArc()) { |
1207 findJoinNode(); |
1305 findJoinNode(); |
1208 bool change = findLeavingArc(); |
1306 bool change = findLeavingArc(); |
1209 changeFlow(change); |
1307 changeFlow(change); |
1210 if (change) { |
1308 if (change) { |
1217 for (int e = _arc_num; e != _arc_num + _node_num; ++e) { |
1315 for (int e = _arc_num; e != _arc_num + _node_num; ++e) { |
1218 if (_flow[e] > 0) return false; |
1316 if (_flow[e] > 0) return false; |
1219 } |
1317 } |
1220 |
1318 |
1221 // Copy flow values to _flow_map |
1319 // Copy flow values to _flow_map |
1222 if (_orig_lower) { |
1320 if (_plower) { |
1223 for (int i = 0; i != _arc_num; ++i) { |
1321 for (int i = 0; i != _arc_num; ++i) { |
1224 Arc e = _arc_ref[i]; |
1322 Arc e = _arc_ref[i]; |
1225 _flow_map->set(e, (*_orig_lower)[e] + _flow[i]); |
1323 _flow_map->set(e, (*_plower)[e] + _flow[i]); |
1226 } |
1324 } |
1227 } else { |
1325 } else { |
1228 for (int i = 0; i != _arc_num; ++i) { |
1326 for (int i = 0; i != _arc_num; ++i) { |
1229 _flow_map->set(_arc_ref[i], _flow[i]); |
1327 _flow_map->set(_arc_ref[i], _flow[i]); |
1230 } |
1328 } |