127 /// The type of the flow map. |
101 /// The type of the flow map. |
128 typedef typename Graph::template EdgeMap<Capacity> FlowMap; |
102 typedef typename Graph::template EdgeMap<Capacity> FlowMap; |
129 /// The type of the potential map. |
103 /// The type of the potential map. |
130 typedef typename Graph::template NodeMap<Cost> PotentialMap; |
104 typedef typename Graph::template NodeMap<Cost> PotentialMap; |
131 |
105 |
132 protected: |
106 public: |
133 |
107 |
|
108 /// Enum type to select the pivot rule used by \ref run(). |
|
109 enum PivotRuleEnum { |
|
110 FIRST_ELIGIBLE_PIVOT, |
|
111 BEST_ELIGIBLE_PIVOT, |
|
112 BLOCK_SEARCH_PIVOT, |
|
113 LIMITED_SEARCH_PIVOT, |
|
114 CANDIDATE_LIST_PIVOT, |
|
115 COMBINED_PIVOT |
|
116 }; |
|
117 |
|
118 private: |
|
119 |
|
120 /// \brief Map adaptor class for handling reduced edge costs. |
|
121 /// |
134 /// Map adaptor class for handling reduced edge costs. |
122 /// Map adaptor class for handling reduced edge costs. |
135 class ReducedCostMap : public MapBase<Edge, Cost> |
123 class ReducedCostMap : public MapBase<Edge, Cost> |
136 { |
124 { |
137 private: |
125 private: |
138 |
126 |
139 const SGraph &gr; |
127 const SGraph &_gr; |
140 const SCostMap &cost_map; |
128 const SCostMap &_cost_map; |
141 const SPotentialMap &pot_map; |
129 const SPotentialMap &_pot_map; |
142 |
130 |
143 public: |
131 public: |
144 |
132 |
145 ReducedCostMap( const SGraph &_gr, |
133 ///\e |
146 const SCostMap &_cm, |
134 ReducedCostMap( const SGraph &gr, |
147 const SPotentialMap &_pm ) : |
135 const SCostMap &cost_map, |
148 gr(_gr), cost_map(_cm), pot_map(_pm) {} |
136 const SPotentialMap &pot_map ) : |
149 |
137 _gr(gr), _cost_map(cost_map), _pot_map(pm) {} |
|
138 |
|
139 ///\e |
150 Cost operator[](const Edge &e) const { |
140 Cost operator[](const Edge &e) const { |
151 return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)]; |
141 return _cost_map[e] + _pot_map[_gr.source(e)] |
|
142 - _pot_map[_gr.target(e)]; |
152 } |
143 } |
153 |
144 |
154 }; //class ReducedCostMap |
145 }; //class ReducedCostMap |
155 |
146 |
156 protected: |
147 private: |
157 |
148 |
158 /// The directed graph the algorithm runs on. |
149 /// \brief Implementation of the "First Eligible" pivot rule for the |
159 SGraph graph; |
150 /// \ref NetworkSimplex "network simplex" algorithm. |
160 /// The original graph. |
151 /// |
161 const Graph &graph_ref; |
152 /// This class implements the "First Eligible" pivot rule |
162 /// The original lower bound map. |
153 /// for the \ref NetworkSimplex "network simplex" algorithm. |
163 const LowerMap *lower; |
154 class FirstEligiblePivotRule |
164 /// The capacity map. |
155 { |
165 SCapacityMap capacity; |
156 private: |
166 /// The cost map. |
157 |
167 SCostMap cost; |
158 NetworkSimplex &_ns; |
168 /// The supply map. |
159 EdgeIt _next_edge; |
169 SSupplyMap supply; |
160 |
170 /// The reduced cost map. |
161 public: |
171 ReducedCostMap red_cost; |
162 |
172 bool valid_supply; |
163 /// Constructor. |
173 |
164 FirstEligiblePivotRule(NetworkSimplex &ns) : |
174 /// The edge map of the current flow. |
165 _ns(ns), _next_edge(ns._graph) {} |
175 SCapacityMap flow; |
166 |
176 /// The potential node map. |
167 /// Finds the next entering edge. |
177 SPotentialMap potential; |
168 bool findEnteringEdge() { |
178 FlowMap flow_result; |
169 for (EdgeIt e = _next_edge; e != INVALID; ++e) { |
179 PotentialMap potential_result; |
170 if (_ns._state[e] * _ns._red_cost[e] < 0) { |
180 |
171 _ns._in_edge = e; |
181 /// Node reference for the original graph. |
172 _next_edge = ++e; |
182 NodeRefMap node_ref; |
173 return true; |
183 /// Edge reference for the original graph. |
174 } |
184 EdgeRefMap edge_ref; |
175 } |
185 |
176 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { |
186 /// The \c depth node map of the spanning tree structure. |
177 if (_ns._state[e] * _ns._red_cost[e] < 0) { |
187 IntNodeMap depth; |
178 _ns._in_edge = e; |
188 /// The \c parent node map of the spanning tree structure. |
179 _next_edge = ++e; |
189 NodeNodeMap parent; |
180 return true; |
190 /// The \c pred_edge node map of the spanning tree structure. |
181 } |
191 EdgeNodeMap pred_edge; |
182 } |
192 /// The \c thread node map of the spanning tree structure. |
183 return false; |
193 NodeNodeMap thread; |
184 } |
194 /// The \c forward node map of the spanning tree structure. |
185 }; //class FirstEligiblePivotRule |
195 BoolNodeMap forward; |
186 |
196 /// The state edge map. |
187 /// \brief Implementation of the "Best Eligible" pivot rule for the |
197 IntEdgeMap state; |
188 /// \ref NetworkSimplex "network simplex" algorithm. |
198 /// The root node of the starting spanning tree. |
189 /// |
199 Node root; |
190 /// This class implements the "Best Eligible" pivot rule |
|
191 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
192 class BestEligiblePivotRule |
|
193 { |
|
194 private: |
|
195 |
|
196 NetworkSimplex &_ns; |
|
197 |
|
198 public: |
|
199 |
|
200 /// Constructor. |
|
201 BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {} |
|
202 |
|
203 /// Finds the next entering edge. |
|
204 bool findEnteringEdge() { |
|
205 Cost min = 0; |
|
206 for (EdgeIt e(_ns._graph); e != INVALID; ++e) { |
|
207 if (_ns._state[e] * _ns._red_cost[e] < min) { |
|
208 min = _ns._state[e] * _ns._red_cost[e]; |
|
209 _ns._in_edge = e; |
|
210 } |
|
211 } |
|
212 return min < 0; |
|
213 } |
|
214 }; //class BestEligiblePivotRule |
|
215 |
|
216 /// \brief Implementation of the "Block Search" pivot rule for the |
|
217 /// \ref NetworkSimplex "network simplex" algorithm. |
|
218 /// |
|
219 /// This class implements the "Block Search" pivot rule |
|
220 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
221 class BlockSearchPivotRule |
|
222 { |
|
223 private: |
|
224 |
|
225 NetworkSimplex &_ns; |
|
226 EdgeIt _next_edge, _min_edge; |
|
227 int _block_size; |
|
228 |
|
229 static const int MIN_BLOCK_SIZE = 10; |
|
230 |
|
231 public: |
|
232 |
|
233 /// Constructor. |
|
234 BlockSearchPivotRule(NetworkSimplex &ns) : |
|
235 _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) |
|
236 { |
|
237 _block_size = 2 * int(sqrt(countEdges(_ns._graph))); |
|
238 if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE; |
|
239 } |
|
240 |
|
241 /// Finds the next entering edge. |
|
242 bool findEnteringEdge() { |
|
243 Cost curr, min = 0; |
|
244 int cnt = 0; |
|
245 for (EdgeIt e = _next_edge; e != INVALID; ++e) { |
|
246 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
|
247 min = curr; |
|
248 _min_edge = e; |
|
249 } |
|
250 if (++cnt == _block_size) { |
|
251 if (min < 0) break; |
|
252 cnt = 0; |
|
253 } |
|
254 } |
|
255 if (min == 0) { |
|
256 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { |
|
257 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
|
258 min = curr; |
|
259 _min_edge = e; |
|
260 } |
|
261 if (++cnt == _block_size) { |
|
262 if (min < 0) break; |
|
263 cnt = 0; |
|
264 } |
|
265 } |
|
266 } |
|
267 _ns._in_edge = _min_edge; |
|
268 _next_edge = ++_min_edge; |
|
269 return min < 0; |
|
270 } |
|
271 }; //class BlockSearchPivotRule |
|
272 |
|
273 /// \brief Implementation of the "Limited Search" pivot rule for the |
|
274 /// \ref NetworkSimplex "network simplex" algorithm. |
|
275 /// |
|
276 /// This class implements the "Limited Search" pivot rule |
|
277 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
278 class LimitedSearchPivotRule |
|
279 { |
|
280 private: |
|
281 |
|
282 NetworkSimplex &_ns; |
|
283 EdgeIt _next_edge, _min_edge; |
|
284 int _sample_size; |
|
285 |
|
286 static const int MIN_SAMPLE_SIZE = 10; |
|
287 static const double SAMPLE_SIZE_FACTOR = 0.0015; |
|
288 |
|
289 public: |
|
290 |
|
291 /// Constructor. |
|
292 LimitedSearchPivotRule(NetworkSimplex &ns) : |
|
293 _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) |
|
294 { |
|
295 _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph)); |
|
296 if (_sample_size < MIN_SAMPLE_SIZE) |
|
297 _sample_size = MIN_SAMPLE_SIZE; |
|
298 } |
|
299 |
|
300 /// Finds the next entering edge. |
|
301 bool findEnteringEdge() { |
|
302 Cost curr, min = 0; |
|
303 int cnt = 0; |
|
304 for (EdgeIt e = _next_edge; e != INVALID; ++e) { |
|
305 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
|
306 min = curr; |
|
307 _min_edge = e; |
|
308 } |
|
309 if (curr < 0 && ++cnt == _sample_size) break; |
|
310 } |
|
311 if (min == 0) { |
|
312 for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { |
|
313 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
|
314 min = curr; |
|
315 _min_edge = e; |
|
316 } |
|
317 if (curr < 0 && ++cnt == _sample_size) break; |
|
318 } |
|
319 } |
|
320 _ns._in_edge = _min_edge; |
|
321 _next_edge = ++_min_edge; |
|
322 return min < 0; |
|
323 } |
|
324 }; //class LimitedSearchPivotRule |
|
325 |
|
326 /// \brief Implementation of the "Candidate List" pivot rule for the |
|
327 /// \ref NetworkSimplex "network simplex" algorithm. |
|
328 /// |
|
329 /// This class implements the "Candidate List" pivot rule |
|
330 /// for the \ref NetworkSimplex "network simplex" algorithm. |
|
331 class CandidateListPivotRule |
|
332 { |
|
333 private: |
|
334 |
|
335 NetworkSimplex &_ns; |
|
336 |
|
337 // The list of candidate edges. |
|
338 std::vector<Edge> _candidates; |
|
339 // The maximum length of the edge list. |
|
340 int _list_length; |
|
341 // The maximum number of minor iterations between two major |
|
342 // itarations. |
|
343 int _minor_limit; |
|
344 |
|
345 int _minor_count; |
|
346 EdgeIt _next_edge; |
|
347 |
|
348 static const double LIST_LENGTH_FACTOR = 0.002; |
|
349 static const double MINOR_LIMIT_FACTOR = 0.1; |
|
350 static const int MIN_LIST_LENGTH = 10; |
|
351 static const int MIN_MINOR_LIMIT = 2; |
|
352 |
|
353 public: |
|
354 |
|
355 /// Constructor. |
|
356 CandidateListPivotRule(NetworkSimplex &ns) : |
|
357 _ns(ns), _next_edge(ns._graph) |
|
358 { |
|
359 int edge_num = countEdges(_ns._graph); |
|
360 _minor_count = 0; |
|
361 _list_length = int(edge_num * LIST_LENGTH_FACTOR); |
|
362 if (_list_length < MIN_LIST_LENGTH) |
|
363 _list_length = MIN_LIST_LENGTH; |
|
364 _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR); |
|
365 if (_minor_limit < MIN_MINOR_LIMIT) |
|
366 _minor_limit = MIN_MINOR_LIMIT; |
|
367 } |
|
368 |
|
369 /// Finds the next entering edge. |
|
370 bool findEnteringEdge() { |
|
371 Cost min, curr; |
|
372 if (_minor_count < _minor_limit && _candidates.size() > 0) { |
|
373 // Minor iteration |
|
374 ++_minor_count; |
|
375 Edge e; |
|
376 min = 0; |
|
377 for (int i = 0; i < int(_candidates.size()); ++i) { |
|
378 e = _candidates[i]; |
|
379 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
|
380 min = curr; |
|
381 _ns._in_edge = e; |
|
382 } |
|
383 } |
|
384 if (min < 0) return true; |
|
385 } |
|
386 |
|
387 // Major iteration |
|
388 _candidates.clear(); |
|
389 EdgeIt e = _next_edge; |
|
390 min = 0; |
|
391 for ( ; e != INVALID; ++e) { |
|
392 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { |
|
393 _candidates.push_back(e); |
|
394 if (curr < min) { |
|
395 min = curr; |
|
396 _ns._in_edge = e; |
|
397 } |
|
398 if (int(_candidates.size()) == _list_length) break; |
|
399 } |
|
400 } |
|
401 if (int(_candidates.size()) < _list_length) { |
|
402 for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) { |
|
403 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { |
|
404 _candidates.push_back(e); |
|
405 if (curr < min) { |
|
406 min = curr; |
|
407 _ns._in_edge = e; |
|
408 } |
|
409 if (int(_candidates.size()) == _list_length) break; |
|
410 } |
|
411 } |
|
412 } |
|
413 if (_candidates.size() == 0) return false; |
|
414 _minor_count = 1; |
|
415 _next_edge = ++e; |
|
416 return true; |
|
417 } |
|
418 }; //class CandidateListPivotRule |
|
419 |
|
420 private: |
|
421 |
|
422 // State constant for edges at their lower bounds. |
|
423 static const int STATE_LOWER = 1; |
|
424 // State constant for edges in the spanning tree. |
|
425 static const int STATE_TREE = 0; |
|
426 // State constant for edges at their upper bounds. |
|
427 static const int STATE_UPPER = -1; |
|
428 |
|
429 // Constant for the combined pivot rule. |
|
430 static const int COMBINED_PIVOT_MAX_DEG = 5; |
|
431 |
|
432 private: |
|
433 |
|
434 // The directed graph the algorithm runs on |
|
435 SGraph _graph; |
|
436 // The original graph |
|
437 const Graph &_graph_ref; |
|
438 // The original lower bound map |
|
439 const LowerMap *_lower; |
|
440 // The capacity map |
|
441 SCapacityMap _capacity; |
|
442 // The cost map |
|
443 SCostMap _cost; |
|
444 // The supply map |
|
445 SSupplyMap _supply; |
|
446 bool _valid_supply; |
|
447 |
|
448 // Edge map of the current flow |
|
449 SCapacityMap _flow; |
|
450 // Node map of the current potentials |
|
451 SPotentialMap _potential; |
|
452 |
|
453 // The depth node map of the spanning tree structure |
|
454 IntNodeMap _depth; |
|
455 // The parent node map of the spanning tree structure |
|
456 NodeNodeMap _parent; |
|
457 // The pred_edge node map of the spanning tree structure |
|
458 EdgeNodeMap _pred_edge; |
|
459 // The thread node map of the spanning tree structure |
|
460 NodeNodeMap _thread; |
|
461 // The forward node map of the spanning tree structure |
|
462 BoolNodeMap _forward; |
|
463 // The state edge map |
|
464 IntEdgeMap _state; |
|
465 // The root node of the starting spanning tree |
|
466 Node _root; |
|
467 |
|
468 // The reduced cost map |
|
469 ReducedCostMap _red_cost; |
|
470 |
|
471 // Members for handling the original graph |
|
472 FlowMap _flow_result; |
|
473 PotentialMap _potential_result; |
|
474 NodeRefMap _node_ref; |
|
475 EdgeRefMap _edge_ref; |
200 |
476 |
201 // The entering edge of the current pivot iteration. |
477 // The entering edge of the current pivot iteration. |
202 Edge in_edge; |
478 Edge _in_edge; |
|
479 |
203 // Temporary nodes used in the current pivot iteration. |
480 // Temporary nodes used in the current pivot iteration. |
204 Node join, u_in, v_in, u_out, v_out; |
481 Node join, u_in, v_in, u_out, v_out; |
205 Node right, first, second, last; |
482 Node right, first, second, last; |
206 Node stem, par_stem, new_stem; |
483 Node stem, par_stem, new_stem; |
207 // The maximum augment amount along the found cycle in the current |
484 // The maximum augment amount along the found cycle in the current |
208 // pivot iteration. |
485 // pivot iteration. |
209 Capacity delta; |
486 Capacity delta; |
210 |
487 |
211 #ifdef EDGE_BLOCK_PIVOT |
|
212 /// The size of blocks for the "Edge Block" pivot rule. |
|
213 int block_size; |
|
214 #endif |
|
215 #ifdef CANDIDATE_LIST_PIVOT |
|
216 /// \brief The list of candidate edges for the "Candidate List" |
|
217 /// pivot rule. |
|
218 std::vector<Edge> candidates; |
|
219 /// \brief The maximum length of the edge list for the |
|
220 /// "Candidate List" pivot rule. |
|
221 int list_length; |
|
222 /// \brief The maximum number of minor iterations between two major |
|
223 /// itarations. |
|
224 int minor_limit; |
|
225 /// \brief The number of minor iterations. |
|
226 int minor_count; |
|
227 #endif |
|
228 #ifdef SORTED_LIST_PIVOT |
|
229 /// \brief The list of candidate edges for the "Sorted List" |
|
230 /// pivot rule. |
|
231 std::vector<Edge> candidates; |
|
232 /// \brief The maximum length of the edge list for the |
|
233 /// "Sorted List" pivot rule. |
|
234 int list_length; |
|
235 int list_index; |
|
236 #endif |
|
237 |
|
238 public : |
488 public : |
239 |
489 |
240 /// \brief General constructor of the class (with lower bounds). |
490 /// \brief General constructor of the class (with lower bounds). |
241 /// |
491 /// |
242 /// General constructor of the class (with lower bounds). |
492 /// General constructor of the class (with lower bounds). |
243 /// |
493 /// |
244 /// \param _graph The directed graph the algorithm runs on. |
494 /// \param graph The directed graph the algorithm runs on. |
245 /// \param _lower The lower bounds of the edges. |
495 /// \param lower The lower bounds of the edges. |
246 /// \param _capacity The capacities (upper bounds) of the edges. |
496 /// \param capacity The capacities (upper bounds) of the edges. |
247 /// \param _cost The cost (length) values of the edges. |
497 /// \param cost The cost (length) values of the edges. |
248 /// \param _supply The supply values of the nodes (signed). |
498 /// \param supply The supply values of the nodes (signed). |
249 NetworkSimplex( const Graph &_graph, |
499 NetworkSimplex( const Graph &graph, |
250 const LowerMap &_lower, |
500 const LowerMap &lower, |
251 const CapacityMap &_capacity, |
501 const CapacityMap &capacity, |
252 const CostMap &_cost, |
502 const CostMap &cost, |
253 const SupplyMap &_supply ) : |
503 const SupplyMap &supply ) : |
254 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), |
504 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), |
255 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
505 _cost(_graph), _supply(_graph), _flow(_graph), |
256 potential_result(_graph), depth(graph), parent(graph), |
506 _potential(_graph), _depth(_graph), _parent(_graph), |
257 pred_edge(graph), thread(graph), forward(graph), state(graph), |
507 _pred_edge(_graph), _thread(_graph), _forward(_graph), |
258 node_ref(graph_ref), edge_ref(graph_ref), |
508 _state(_graph), _red_cost(_graph, _cost, _potential), |
259 red_cost(graph, cost, potential) |
509 _flow_result(graph), _potential_result(graph), |
|
510 _node_ref(graph), _edge_ref(graph) |
260 { |
511 { |
261 // Checking the sum of supply values |
512 // Checking the sum of supply values |
262 Supply sum = 0; |
513 Supply sum = 0; |
263 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
514 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
264 sum += _supply[n]; |
515 sum += supply[n]; |
265 if (!(valid_supply = sum == 0)) return; |
516 if (!(_valid_supply = sum == 0)) return; |
266 |
517 |
267 // Copying graph_ref to graph |
518 // Copying _graph_ref to _graph |
268 graph.reserveNode(countNodes(graph_ref) + 1); |
519 _graph.reserveNode(countNodes(_graph_ref) + 1); |
269 graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref)); |
520 _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); |
270 copyGraph(graph, graph_ref) |
521 copyGraph(_graph, _graph_ref) |
271 .edgeMap(cost, _cost) |
522 .edgeMap(_cost, cost) |
272 .nodeRef(node_ref) |
523 .nodeRef(_node_ref) |
273 .edgeRef(edge_ref) |
524 .edgeRef(_edge_ref) |
274 .run(); |
525 .run(); |
275 |
526 |
276 // Removing non-zero lower bounds |
527 // Removing non-zero lower bounds |
277 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { |
528 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { |
278 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; |
529 _capacity[_edge_ref[e]] = capacity[e] - lower[e]; |
279 } |
530 } |
280 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { |
531 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { |
281 Supply s = _supply[n]; |
532 Supply s = supply[n]; |
282 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) |
533 for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
283 s += _lower[e]; |
534 s += lower[e]; |
284 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) |
535 for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
285 s -= _lower[e]; |
536 s -= lower[e]; |
286 supply[node_ref[n]] = s; |
537 _supply[_node_ref[n]] = s; |
287 } |
538 } |
288 } |
539 } |
289 |
540 |
290 /// \brief General constructor of the class (without lower bounds). |
541 /// \brief General constructor of the class (without lower bounds). |
291 /// |
542 /// |
292 /// General constructor of the class (without lower bounds). |
543 /// General constructor of the class (without lower bounds). |
293 /// |
544 /// |
294 /// \param _graph The directed graph the algorithm runs on. |
545 /// \param graph The directed graph the algorithm runs on. |
295 /// \param _capacity The capacities (upper bounds) of the edges. |
546 /// \param capacity The capacities (upper bounds) of the edges. |
296 /// \param _cost The cost (length) values of the edges. |
547 /// \param cost The cost (length) values of the edges. |
297 /// \param _supply The supply values of the nodes (signed). |
548 /// \param supply The supply values of the nodes (signed). |
298 NetworkSimplex( const Graph &_graph, |
549 NetworkSimplex( const Graph &graph, |
299 const CapacityMap &_capacity, |
550 const CapacityMap &capacity, |
300 const CostMap &_cost, |
551 const CostMap &cost, |
301 const SupplyMap &_supply ) : |
552 const SupplyMap &supply ) : |
302 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), |
553 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), |
303 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
554 _cost(_graph), _supply(_graph), _flow(_graph), |
304 potential_result(_graph), depth(graph), parent(graph), |
555 _potential(_graph), _depth(_graph), _parent(_graph), |
305 pred_edge(graph), thread(graph), forward(graph), state(graph), |
556 _pred_edge(_graph), _thread(_graph), _forward(_graph), |
306 node_ref(graph_ref), edge_ref(graph_ref), |
557 _state(_graph), _red_cost(_graph, _cost, _potential), |
307 red_cost(graph, cost, potential) |
558 _flow_result(graph), _potential_result(graph), |
|
559 _node_ref(graph), _edge_ref(graph) |
308 { |
560 { |
309 // Checking the sum of supply values |
561 // Checking the sum of supply values |
310 Supply sum = 0; |
562 Supply sum = 0; |
311 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
563 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
312 sum += _supply[n]; |
564 sum += supply[n]; |
313 if (!(valid_supply = sum == 0)) return; |
565 if (!(_valid_supply = sum == 0)) return; |
314 |
566 |
315 // Copying graph_ref to graph |
567 // Copying _graph_ref to graph |
316 copyGraph(graph, graph_ref) |
568 copyGraph(_graph, _graph_ref) |
317 .edgeMap(capacity, _capacity) |
569 .edgeMap(_capacity, capacity) |
318 .edgeMap(cost, _cost) |
570 .edgeMap(_cost, cost) |
319 .nodeMap(supply, _supply) |
571 .nodeMap(_supply, supply) |
320 .nodeRef(node_ref) |
572 .nodeRef(_node_ref) |
321 .edgeRef(edge_ref) |
573 .edgeRef(_edge_ref) |
322 .run(); |
574 .run(); |
323 } |
575 } |
324 |
576 |
325 /// \brief Simple constructor of the class (with lower bounds). |
577 /// \brief Simple constructor of the class (with lower bounds). |
326 /// |
578 /// |
327 /// Simple constructor of the class (with lower bounds). |
579 /// Simple constructor of the class (with lower bounds). |
328 /// |
580 /// |
329 /// \param _graph The directed graph the algorithm runs on. |
581 /// \param graph The directed graph the algorithm runs on. |
330 /// \param _lower The lower bounds of the edges. |
582 /// \param lower The lower bounds of the edges. |
331 /// \param _capacity The capacities (upper bounds) of the edges. |
583 /// \param capacity The capacities (upper bounds) of the edges. |
332 /// \param _cost The cost (length) values of the edges. |
584 /// \param cost The cost (length) values of the edges. |
333 /// \param _s The source node. |
585 /// \param s The source node. |
334 /// \param _t The target node. |
586 /// \param t The target node. |
335 /// \param _flow_value The required amount of flow from node \c _s |
587 /// \param flow_value The required amount of flow from node \c s |
336 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). |
588 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
337 NetworkSimplex( const Graph &_graph, |
589 NetworkSimplex( const Graph &graph, |
338 const LowerMap &_lower, |
590 const LowerMap &lower, |
339 const CapacityMap &_capacity, |
591 const CapacityMap &capacity, |
340 const CostMap &_cost, |
592 const CostMap &cost, |
341 typename Graph::Node _s, |
593 typename Graph::Node s, |
342 typename Graph::Node _t, |
594 typename Graph::Node t, |
343 typename SupplyMap::Value _flow_value ) : |
595 typename SupplyMap::Value flow_value ) : |
344 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), |
596 _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), |
345 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
597 _cost(_graph), _supply(_graph), _flow(_graph), |
346 potential_result(_graph), depth(graph), parent(graph), |
598 _potential(_graph), _depth(_graph), _parent(_graph), |
347 pred_edge(graph), thread(graph), forward(graph), state(graph), |
599 _pred_edge(_graph), _thread(_graph), _forward(_graph), |
348 node_ref(graph_ref), edge_ref(graph_ref), |
600 _state(_graph), _red_cost(_graph, _cost, _potential), |
349 red_cost(graph, cost, potential) |
601 _flow_result(graph), _potential_result(graph), |
|
602 _node_ref(graph), _edge_ref(graph) |
350 { |
603 { |
351 // Copying graph_ref to graph |
604 // Copying _graph_ref to graph |
352 copyGraph(graph, graph_ref) |
605 copyGraph(_graph, _graph_ref) |
353 .edgeMap(cost, _cost) |
606 .edgeMap(_cost, cost) |
354 .nodeRef(node_ref) |
607 .nodeRef(_node_ref) |
355 .edgeRef(edge_ref) |
608 .edgeRef(_edge_ref) |
356 .run(); |
609 .run(); |
357 |
610 |
358 // Removing non-zero lower bounds |
611 // Removing non-zero lower bounds |
359 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { |
612 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { |
360 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; |
613 _capacity[_edge_ref[e]] = capacity[e] - lower[e]; |
361 } |
614 } |
362 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { |
615 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { |
363 Supply s = 0; |
616 Supply sum = 0; |
364 if (n == _s) s = _flow_value; |
617 if (n == s) sum = flow_value; |
365 if (n == _t) s = -_flow_value; |
618 if (n == t) sum = -flow_value; |
366 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) |
619 for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
367 s += _lower[e]; |
620 sum += lower[e]; |
368 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) |
621 for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
369 s -= _lower[e]; |
622 sum -= lower[e]; |
370 supply[node_ref[n]] = s; |
623 _supply[_node_ref[n]] = sum; |
371 } |
624 } |
372 valid_supply = true; |
625 _valid_supply = true; |
373 } |
626 } |
374 |
627 |
375 /// \brief Simple constructor of the class (without lower bounds). |
628 /// \brief Simple constructor of the class (without lower bounds). |
376 /// |
629 /// |
377 /// Simple constructor of the class (without lower bounds). |
630 /// Simple constructor of the class (without lower bounds). |
378 /// |
631 /// |
379 /// \param _graph The directed graph the algorithm runs on. |
632 /// \param graph The directed graph the algorithm runs on. |
380 /// \param _capacity The capacities (upper bounds) of the edges. |
633 /// \param capacity The capacities (upper bounds) of the edges. |
381 /// \param _cost The cost (length) values of the edges. |
634 /// \param cost The cost (length) values of the edges. |
382 /// \param _s The source node. |
635 /// \param s The source node. |
383 /// \param _t The target node. |
636 /// \param t The target node. |
384 /// \param _flow_value The required amount of flow from node \c _s |
637 /// \param flow_value The required amount of flow from node \c s |
385 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). |
638 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
386 NetworkSimplex( const Graph &_graph, |
639 NetworkSimplex( const Graph &graph, |
387 const CapacityMap &_capacity, |
640 const CapacityMap &capacity, |
388 const CostMap &_cost, |
641 const CostMap &cost, |
389 typename Graph::Node _s, |
642 typename Graph::Node s, |
390 typename Graph::Node _t, |
643 typename Graph::Node t, |
391 typename SupplyMap::Value _flow_value ) : |
644 typename SupplyMap::Value flow_value ) : |
392 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), |
645 _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), |
393 supply(graph, 0), flow(graph), flow_result(_graph), potential(graph), |
646 _cost(_graph), _supply(_graph, 0), _flow(_graph), |
394 potential_result(_graph), depth(graph), parent(graph), |
647 _potential(_graph), _depth(_graph), _parent(_graph), |
395 pred_edge(graph), thread(graph), forward(graph), state(graph), |
648 _pred_edge(_graph), _thread(_graph), _forward(_graph), |
396 node_ref(graph_ref), edge_ref(graph_ref), |
649 _state(_graph), _red_cost(_graph, _cost, _potential), |
397 red_cost(graph, cost, potential) |
650 _flow_result(graph), _potential_result(graph), |
|
651 _node_ref(graph), _edge_ref(graph) |
398 { |
652 { |
399 // Copying graph_ref to graph |
653 // Copying _graph_ref to graph |
400 copyGraph(graph, graph_ref) |
654 copyGraph(_graph, _graph_ref) |
401 .edgeMap(capacity, _capacity) |
655 .edgeMap(_capacity, capacity) |
402 .edgeMap(cost, _cost) |
656 .edgeMap(_cost, cost) |
403 .nodeRef(node_ref) |
657 .nodeRef(_node_ref) |
404 .edgeRef(edge_ref) |
658 .edgeRef(_edge_ref) |
405 .run(); |
659 .run(); |
406 supply[node_ref[_s]] = _flow_value; |
660 _supply[_node_ref[s]] = flow_value; |
407 supply[node_ref[_t]] = -_flow_value; |
661 _supply[_node_ref[t]] = -flow_value; |
408 valid_supply = true; |
662 _valid_supply = true; |
409 } |
663 } |
410 |
664 |
411 /// \brief Runs the algorithm. |
665 /// \brief Runs the algorithm. |
412 /// |
666 /// |
413 /// Runs the algorithm. |
667 /// Runs the algorithm. |
414 /// |
668 /// |
|
669 /// \param pivot_rule The pivot rule that is used during the |
|
670 /// algorithm. |
|
671 /// |
|
672 /// The available pivot rules: |
|
673 /// |
|
674 /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in |
|
675 /// a wraparound fashion in every iteration |
|
676 /// (\ref FirstEligiblePivotRule). |
|
677 /// |
|
678 /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in |
|
679 /// every iteration (\ref BestEligiblePivotRule). |
|
680 /// |
|
681 /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in |
|
682 /// every iteration in a wraparound fashion and the best eligible |
|
683 /// edge is selected from this block (\ref BlockSearchPivotRule). |
|
684 /// |
|
685 /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are |
|
686 /// examined in every iteration in a wraparound fashion and the best |
|
687 /// one is selected from them (\ref LimitedSearchPivotRule). |
|
688 /// |
|
689 /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is |
|
690 /// built from eligible edges and it is used for edge selection in |
|
691 /// the following minor iterations (\ref CandidateListPivotRule). |
|
692 /// |
|
693 /// - COMBINED_PIVOT This is a combined version of the two fastest |
|
694 /// pivot rules. |
|
695 /// For rather sparse graphs \ref LimitedSearchPivotRule |
|
696 /// "Limited Search" implementation is used, otherwise |
|
697 /// \ref BlockSearchPivotRule "Block Search" pivot rule is used. |
|
698 /// According to our benchmark tests this combined method is the |
|
699 /// most efficient. |
|
700 /// |
415 /// \return \c true if a feasible flow can be found. |
701 /// \return \c true if a feasible flow can be found. |
416 bool run() { |
702 bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) { |
417 return init() && start(); |
703 return init() && start(pivot_rule); |
418 } |
704 } |
419 |
705 |
420 /// \brief Returns a const reference to the flow map. |
706 /// \brief Returns a const reference to the edge map storing the |
421 /// |
707 /// found flow. |
422 /// Returns a const reference to the flow map. |
708 /// |
|
709 /// Returns a const reference to the edge map storing the found flow. |
423 /// |
710 /// |
424 /// \pre \ref run() must be called before using this function. |
711 /// \pre \ref run() must be called before using this function. |
425 const FlowMap& flowMap() const { |
712 const FlowMap& flowMap() const { |
426 return flow_result; |
713 return _flow_result; |
427 } |
714 } |
428 |
715 |
429 /// \brief Returns a const reference to the potential map (the dual |
716 /// \brief Returns a const reference to the node map storing the |
430 /// solution). |
717 /// found potentials (the dual solution). |
431 /// |
718 /// |
432 /// Returns a const reference to the potential map (the dual |
719 /// Returns a const reference to the node map storing the found |
433 /// solution). |
720 /// potentials (the dual solution). |
434 /// |
721 /// |
435 /// \pre \ref run() must be called before using this function. |
722 /// \pre \ref run() must be called before using this function. |
436 const PotentialMap& potentialMap() const { |
723 const PotentialMap& potentialMap() const { |
437 return potential_result; |
724 return _potential_result; |
438 } |
725 } |
439 |
726 |
440 /// \brief Returns the total cost of the found flow. |
727 /// \brief Returns the total cost of the found flow. |
441 /// |
728 /// |
442 /// Returns the total cost of the found flow. The complexity of the |
729 /// Returns the total cost of the found flow. The complexity of the |
443 /// function is \f$ O(e) \f$. |
730 /// function is \f$ O(e) \f$. |
444 /// |
731 /// |
445 /// \pre \ref run() must be called before using this function. |
732 /// \pre \ref run() must be called before using this function. |
446 Cost totalCost() const { |
733 Cost totalCost() const { |
447 Cost c = 0; |
734 Cost c = 0; |
448 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) |
735 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
449 c += flow_result[e] * cost[edge_ref[e]]; |
736 c += _flow_result[e] * _cost[_edge_ref[e]]; |
450 return c; |
737 return c; |
451 } |
738 } |
452 |
739 |
453 protected: |
740 private: |
454 |
741 |
455 /// \brief Extends the underlaying graph and initializes all the |
742 /// \brief Extends the underlaying graph and initializes all the |
456 /// node and edge maps. |
743 /// node and edge maps. |
457 bool init() { |
744 bool init() { |
458 if (!valid_supply) return false; |
745 if (!_valid_supply) return false; |
459 |
746 |
460 // Initializing state and flow maps |
747 // Initializing state and flow maps |
461 for (EdgeIt e(graph); e != INVALID; ++e) { |
748 for (EdgeIt e(_graph); e != INVALID; ++e) { |
462 flow[e] = 0; |
749 _flow[e] = 0; |
463 state[e] = LOWER; |
750 _state[e] = STATE_LOWER; |
464 } |
751 } |
465 |
752 |
466 // Adding an artificial root node to the graph |
753 // Adding an artificial root node to the graph |
467 root = graph.addNode(); |
754 _root = _graph.addNode(); |
468 parent[root] = INVALID; |
755 _parent[_root] = INVALID; |
469 pred_edge[root] = INVALID; |
756 _pred_edge[_root] = INVALID; |
470 depth[root] = 0; |
757 _depth[_root] = 0; |
471 supply[root] = 0; |
758 _supply[_root] = 0; |
472 potential[root] = 0; |
759 _potential[_root] = 0; |
473 |
760 |
474 // Adding artificial edges to the graph and initializing the node |
761 // Adding artificial edges to the graph and initializing the node |
475 // maps of the spanning tree data structure |
762 // maps of the spanning tree data structure |
476 Supply sum = 0; |
763 Node last = _root; |
477 Node last = root; |
|
478 Edge e; |
764 Edge e; |
479 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
765 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
480 for (NodeIt u(graph); u != INVALID; ++u) { |
766 for (NodeIt u(_graph); u != INVALID; ++u) { |
481 if (u == root) continue; |
767 if (u == _root) continue; |
482 thread[last] = u; |
768 _thread[last] = u; |
483 last = u; |
769 last = u; |
484 parent[u] = root; |
770 _parent[u] = _root; |
485 depth[u] = 1; |
771 _depth[u] = 1; |
486 sum += supply[u]; |
772 if (_supply[u] >= 0) { |
487 if (supply[u] >= 0) { |
773 e = _graph.addEdge(u, _root); |
488 e = graph.addEdge(u, root); |
774 _flow[e] = _supply[u]; |
489 flow[e] = supply[u]; |
775 _forward[u] = true; |
490 forward[u] = true; |
776 _potential[u] = -max_cost; |
491 potential[u] = max_cost; |
|
492 } else { |
777 } else { |
493 e = graph.addEdge(root, u); |
778 e = _graph.addEdge(_root, u); |
494 flow[e] = -supply[u]; |
779 _flow[e] = -_supply[u]; |
495 forward[u] = false; |
780 _forward[u] = false; |
496 potential[u] = -max_cost; |
781 _potential[u] = max_cost; |
497 } |
782 } |
498 cost[e] = max_cost; |
783 _cost[e] = max_cost; |
499 capacity[e] = std::numeric_limits<Capacity>::max(); |
784 _capacity[e] = std::numeric_limits<Capacity>::max(); |
500 state[e] = TREE; |
785 _state[e] = STATE_TREE; |
501 pred_edge[u] = e; |
786 _pred_edge[u] = e; |
502 } |
787 } |
503 thread[last] = root; |
788 _thread[last] = _root; |
504 |
789 |
505 #ifdef EDGE_BLOCK_PIVOT |
|
506 // Initializing block_size for the edge block pivot rule |
|
507 int edge_num = countEdges(graph); |
|
508 block_size = 2 * int(sqrt(countEdges(graph))); |
|
509 if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE; |
|
510 #endif |
|
511 #ifdef CANDIDATE_LIST_PIVOT |
|
512 int edge_num = countEdges(graph); |
|
513 minor_count = 0; |
|
514 list_length = edge_num / LIST_LENGTH_DIV; |
|
515 minor_limit = edge_num / MINOR_LIMIT_DIV; |
|
516 #endif |
|
517 #ifdef SORTED_LIST_PIVOT |
|
518 int edge_num = countEdges(graph); |
|
519 list_index = 0; |
|
520 list_length = edge_num / LIST_LENGTH_DIV; |
|
521 #endif |
|
522 |
|
523 return sum == 0; |
|
524 } |
|
525 |
|
526 #ifdef FIRST_ELIGIBLE_PIVOT |
|
527 /// \brief Finds entering edge according to the "First Eligible" |
|
528 /// pivot rule. |
|
529 bool findEnteringEdge(EdgeIt &next_edge) { |
|
530 for (EdgeIt e = next_edge; e != INVALID; ++e) { |
|
531 if (state[e] * red_cost[e] < 0) { |
|
532 in_edge = e; |
|
533 next_edge = ++e; |
|
534 return true; |
|
535 } |
|
536 } |
|
537 for (EdgeIt e(graph); e != next_edge; ++e) { |
|
538 if (state[e] * red_cost[e] < 0) { |
|
539 in_edge = e; |
|
540 next_edge = ++e; |
|
541 return true; |
|
542 } |
|
543 } |
|
544 return false; |
|
545 } |
|
546 #endif |
|
547 |
|
548 #ifdef BEST_ELIGIBLE_PIVOT |
|
549 /// \brief Finds entering edge according to the "Best Eligible" |
|
550 /// pivot rule. |
|
551 bool findEnteringEdge() { |
|
552 Cost min = 0; |
|
553 for (EdgeIt e(graph); e != INVALID; ++e) { |
|
554 if (state[e] * red_cost[e] < min) { |
|
555 min = state[e] * red_cost[e]; |
|
556 in_edge = e; |
|
557 } |
|
558 } |
|
559 return min < 0; |
|
560 } |
|
561 #endif |
|
562 |
|
563 #ifdef EDGE_BLOCK_PIVOT |
|
564 /// \brief Finds entering edge according to the "Edge Block" |
|
565 /// pivot rule. |
|
566 bool findEnteringEdge(EdgeIt &next_edge) { |
|
567 // Performing edge block selection |
|
568 Cost curr, min = 0; |
|
569 EdgeIt min_edge(graph); |
|
570 int cnt = 0; |
|
571 for (EdgeIt e = next_edge; e != INVALID; ++e) { |
|
572 if ((curr = state[e] * red_cost[e]) < min) { |
|
573 min = curr; |
|
574 min_edge = e; |
|
575 } |
|
576 if (++cnt == block_size) { |
|
577 if (min < 0) break; |
|
578 cnt = 0; |
|
579 } |
|
580 } |
|
581 if (!(min < 0)) { |
|
582 for (EdgeIt e(graph); e != next_edge; ++e) { |
|
583 if ((curr = state[e] * red_cost[e]) < min) { |
|
584 min = curr; |
|
585 min_edge = e; |
|
586 } |
|
587 if (++cnt == block_size) { |
|
588 if (min < 0) break; |
|
589 cnt = 0; |
|
590 } |
|
591 } |
|
592 } |
|
593 in_edge = min_edge; |
|
594 if ((next_edge = ++min_edge) == INVALID) |
|
595 next_edge = EdgeIt(graph); |
|
596 return min < 0; |
|
597 } |
|
598 #endif |
|
599 |
|
600 #ifdef CANDIDATE_LIST_PIVOT |
|
601 /// \brief Finds entering edge according to the "Candidate List" |
|
602 /// pivot rule. |
|
603 bool findEnteringEdge() { |
|
604 typedef typename std::vector<Edge>::iterator ListIt; |
|
605 |
|
606 if (minor_count >= minor_limit || candidates.size() == 0) { |
|
607 // Major iteration |
|
608 candidates.clear(); |
|
609 for (EdgeIt e(graph); e != INVALID; ++e) { |
|
610 if (state[e] * red_cost[e] < 0) { |
|
611 candidates.push_back(e); |
|
612 if (candidates.size() == list_length) break; |
|
613 } |
|
614 } |
|
615 if (candidates.size() == 0) return false; |
|
616 } |
|
617 |
|
618 // Minor iteration |
|
619 ++minor_count; |
|
620 Cost min = 0; |
|
621 Edge e; |
|
622 for (int i = 0; i < candidates.size(); ++i) { |
|
623 e = candidates[i]; |
|
624 if (state[e] * red_cost[e] < min) { |
|
625 min = state[e] * red_cost[e]; |
|
626 in_edge = e; |
|
627 } |
|
628 } |
|
629 return true; |
790 return true; |
630 } |
791 } |
631 #endif |
792 |
632 |
793 /// Finds the join node. |
633 #ifdef SORTED_LIST_PIVOT |
|
634 /// \brief Functor class to compare edges during sort of the |
|
635 /// candidate list. |
|
636 class SortFunc |
|
637 { |
|
638 private: |
|
639 const IntEdgeMap &st; |
|
640 const ReducedCostMap &rc; |
|
641 public: |
|
642 SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) : |
|
643 st(_st), rc(_rc) {} |
|
644 bool operator()(const Edge &e1, const Edge &e2) { |
|
645 return st[e1] * rc[e1] < st[e2] * rc[e2]; |
|
646 } |
|
647 }; |
|
648 |
|
649 /// \brief Finds entering edge according to the "Sorted List" |
|
650 /// pivot rule. |
|
651 bool findEnteringEdge() { |
|
652 static SortFunc sort_func(state, red_cost); |
|
653 |
|
654 // Minor iteration |
|
655 while (list_index < candidates.size()) { |
|
656 in_edge = candidates[list_index++]; |
|
657 if (state[in_edge] * red_cost[in_edge] < 0) return true; |
|
658 } |
|
659 |
|
660 // Major iteration |
|
661 candidates.clear(); |
|
662 Cost curr, min = 0; |
|
663 for (EdgeIt e(graph); e != INVALID; ++e) { |
|
664 if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) { |
|
665 candidates.push_back(e); |
|
666 if (curr < min) min = curr; |
|
667 if (candidates.size() == list_length) break; |
|
668 } |
|
669 } |
|
670 if (candidates.size() == 0) return false; |
|
671 sort(candidates.begin(), candidates.end(), sort_func); |
|
672 in_edge = candidates[0]; |
|
673 list_index = 1; |
|
674 return true; |
|
675 } |
|
676 #endif |
|
677 |
|
678 /// \brief Finds the join node. |
|
679 Node findJoinNode() { |
794 Node findJoinNode() { |
680 // Finding the join node |
795 Node u = _graph.source(_in_edge); |
681 Node u = graph.source(in_edge); |
796 Node v = _graph.target(_in_edge); |
682 Node v = graph.target(in_edge); |
|
683 while (u != v) { |
797 while (u != v) { |
684 if (depth[u] == depth[v]) { |
798 if (_depth[u] == _depth[v]) { |
685 u = parent[u]; |
799 u = _parent[u]; |
686 v = parent[v]; |
800 v = _parent[v]; |
687 } |
801 } |
688 else if (depth[u] > depth[v]) u = parent[u]; |
802 else if (_depth[u] > _depth[v]) u = _parent[u]; |
689 else v = parent[v]; |
803 else v = _parent[v]; |
690 } |
804 } |
691 return u; |
805 return u; |
692 } |
806 } |
693 |
807 |
694 /// \brief Finds the leaving edge of the cycle. Returns \c true if |
808 /// \brief Finds the leaving edge of the cycle. Returns \c true if |
695 /// the leaving edge is not the same as the entering edge. |
809 /// the leaving edge is not the same as the entering edge. |
696 bool findLeavingEdge() { |
810 bool findLeavingEdge() { |
697 // Initializing first and second nodes according to the direction |
811 // Initializing first and second nodes according to the direction |
698 // of the cycle |
812 // of the cycle |
699 if (state[in_edge] == LOWER) { |
813 if (_state[_in_edge] == STATE_LOWER) { |
700 first = graph.source(in_edge); |
814 first = _graph.source(_in_edge); |
701 second = graph.target(in_edge); |
815 second = _graph.target(_in_edge); |
702 } else { |
816 } else { |
703 first = graph.target(in_edge); |
817 first = _graph.target(_in_edge); |
704 second = graph.source(in_edge); |
818 second = _graph.source(_in_edge); |
705 } |
819 } |
706 delta = capacity[in_edge]; |
820 delta = _capacity[_in_edge]; |
707 bool result = false; |
821 bool result = false; |
708 Capacity d; |
822 Capacity d; |
709 Edge e; |
823 Edge e; |
710 |
824 |
711 // Searching the cycle along the path form the first node to the |
825 // Searching the cycle along the path form the first node to the |
712 // root node |
826 // root node |
713 for (Node u = first; u != join; u = parent[u]) { |
827 for (Node u = first; u != join; u = _parent[u]) { |
714 e = pred_edge[u]; |
828 e = _pred_edge[u]; |
715 d = forward[u] ? flow[e] : capacity[e] - flow[e]; |
829 d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e]; |
716 if (d < delta) { |
830 if (d < delta) { |
717 delta = d; |
831 delta = d; |
718 u_out = u; |
832 u_out = u; |
719 u_in = first; |
833 u_in = first; |
720 v_in = second; |
834 v_in = second; |
721 result = true; |
835 result = true; |
722 } |
836 } |
723 } |
837 } |
724 // Searching the cycle along the path form the second node to the |
838 // Searching the cycle along the path form the second node to the |
725 // root node |
839 // root node |
726 for (Node u = second; u != join; u = parent[u]) { |
840 for (Node u = second; u != join; u = _parent[u]) { |
727 e = pred_edge[u]; |
841 e = _pred_edge[u]; |
728 d = forward[u] ? capacity[e] - flow[e] : flow[e]; |
842 d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e]; |
729 if (d <= delta) { |
843 if (d <= delta) { |
730 delta = d; |
844 delta = d; |
731 u_out = u; |
845 u_out = u; |
732 u_in = second; |
846 u_in = second; |
733 v_in = first; |
847 v_in = first; |
735 } |
849 } |
736 } |
850 } |
737 return result; |
851 return result; |
738 } |
852 } |
739 |
853 |
740 /// \brief Changes \c flow and \c state edge maps. |
854 /// Changes \c flow and \c state edge maps. |
741 void changeFlows(bool change) { |
855 void changeFlows(bool change) { |
742 // Augmenting along the cycle |
856 // Augmenting along the cycle |
743 if (delta > 0) { |
857 if (delta > 0) { |
744 Capacity val = state[in_edge] * delta; |
858 Capacity val = _state[_in_edge] * delta; |
745 flow[in_edge] += val; |
859 _flow[_in_edge] += val; |
746 for (Node u = graph.source(in_edge); u != join; u = parent[u]) { |
860 for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) { |
747 flow[pred_edge[u]] += forward[u] ? -val : val; |
861 _flow[_pred_edge[u]] += _forward[u] ? -val : val; |
748 } |
862 } |
749 for (Node u = graph.target(in_edge); u != join; u = parent[u]) { |
863 for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) { |
750 flow[pred_edge[u]] += forward[u] ? val : -val; |
864 _flow[_pred_edge[u]] += _forward[u] ? val : -val; |
751 } |
865 } |
752 } |
866 } |
753 // Updating the state of the entering and leaving edges |
867 // Updating the state of the entering and leaving edges |
754 if (change) { |
868 if (change) { |
755 state[in_edge] = TREE; |
869 _state[_in_edge] = STATE_TREE; |
756 state[pred_edge[u_out]] = |
870 _state[_pred_edge[u_out]] = |
757 (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER; |
871 (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; |
758 } else { |
872 } else { |
759 state[in_edge] = -state[in_edge]; |
873 _state[_in_edge] = -_state[_in_edge]; |
760 } |
874 } |
761 } |
875 } |
762 |
876 |
763 /// \brief Updates \c thread and \c parent node maps. |
877 /// Updates \c thread and \c parent node maps. |
764 void updateThreadParent() { |
878 void updateThreadParent() { |
765 Node u; |
879 Node u; |
766 v_out = parent[u_out]; |
880 v_out = _parent[u_out]; |
767 |
881 |
768 // Handling the case when join and v_out coincide |
882 // Handling the case when join and v_out coincide |
769 bool par_first = false; |
883 bool par_first = false; |
770 if (join == v_out) { |
884 if (join == v_out) { |
771 for (u = join; u != u_in && u != v_in; u = thread[u]) ; |
885 for (u = join; u != u_in && u != v_in; u = _thread[u]) ; |
772 if (u == v_in) { |
886 if (u == v_in) { |
773 par_first = true; |
887 par_first = true; |
774 while (thread[u] != u_out) u = thread[u]; |
888 while (_thread[u] != u_out) u = _thread[u]; |
775 first = u; |
889 first = u; |
776 } |
890 } |
777 } |
891 } |
778 |
892 |
779 // Finding the last successor of u_in (u) and the node after it |
893 // Finding the last successor of u_in (u) and the node after it |
780 // (right) according to the thread index |
894 // (right) according to the thread index |
781 for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ; |
895 for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ; |
782 right = thread[u]; |
896 right = _thread[u]; |
783 if (thread[v_in] == u_out) { |
897 if (_thread[v_in] == u_out) { |
784 for (last = u; depth[last] > depth[u_out]; last = thread[last]) ; |
898 for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ; |
785 if (last == u_out) last = thread[last]; |
899 if (last == u_out) last = _thread[last]; |
786 } |
900 } |
787 else last = thread[v_in]; |
901 else last = _thread[v_in]; |
788 |
902 |
789 // Updating stem nodes |
903 // Updating stem nodes |
790 thread[v_in] = stem = u_in; |
904 _thread[v_in] = stem = u_in; |
791 par_stem = v_in; |
905 par_stem = v_in; |
792 while (stem != u_out) { |
906 while (stem != u_out) { |
793 thread[u] = new_stem = parent[stem]; |
907 _thread[u] = new_stem = _parent[stem]; |
794 |
908 |
795 // Finding the node just before the stem node (u) according to |
909 // Finding the node just before the stem node (u) according to |
796 // the original thread index |
910 // the original thread index |
797 for (u = new_stem; thread[u] != stem; u = thread[u]) ; |
911 for (u = new_stem; _thread[u] != stem; u = _thread[u]) ; |
798 thread[u] = right; |
912 _thread[u] = right; |
799 |
913 |
800 // Changing the parent node of stem and shifting stem and |
914 // Changing the parent node of stem and shifting stem and |
801 // par_stem nodes |
915 // par_stem nodes |
802 parent[stem] = par_stem; |
916 _parent[stem] = par_stem; |
803 par_stem = stem; |
917 par_stem = stem; |
804 stem = new_stem; |
918 stem = new_stem; |
805 |
919 |
806 // Finding the last successor of stem (u) and the node after it |
920 // Finding the last successor of stem (u) and the node after it |
807 // (right) according to the thread index |
921 // (right) according to the thread index |
808 for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ; |
922 for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ; |
809 right = thread[u]; |
923 right = _thread[u]; |
810 } |
924 } |
811 parent[u_out] = par_stem; |
925 _parent[u_out] = par_stem; |
812 thread[u] = last; |
926 _thread[u] = last; |
813 |
927 |
814 if (join == v_out && par_first) { |
928 if (join == v_out && par_first) { |
815 if (first != v_in) thread[first] = right; |
929 if (first != v_in) _thread[first] = right; |
816 } else { |
930 } else { |
817 for (u = v_out; thread[u] != u_out; u = thread[u]) ; |
931 for (u = v_out; _thread[u] != u_out; u = _thread[u]) ; |
818 thread[u] = right; |
932 _thread[u] = right; |
819 } |
933 } |
820 } |
934 } |
821 |
935 |
822 /// \brief Updates \c pred_edge and \c forward node maps. |
936 /// Updates \c pred_edge and \c forward node maps. |
823 void updatePredEdge() { |
937 void updatePredEdge() { |
824 Node u = u_out, v; |
938 Node u = u_out, v; |
825 while (u != u_in) { |
939 while (u != u_in) { |
826 v = parent[u]; |
940 v = _parent[u]; |
827 pred_edge[u] = pred_edge[v]; |
941 _pred_edge[u] = _pred_edge[v]; |
828 forward[u] = !forward[v]; |
942 _forward[u] = !_forward[v]; |
829 u = v; |
943 u = v; |
830 } |
944 } |
831 pred_edge[u_in] = in_edge; |
945 _pred_edge[u_in] = _in_edge; |
832 forward[u_in] = (u_in == graph.source(in_edge)); |
946 _forward[u_in] = (u_in == _graph.source(_in_edge)); |
833 } |
947 } |
834 |
948 |
835 /// \brief Updates \c depth and \c potential node maps. |
949 /// Updates \c depth and \c potential node maps. |
836 void updateDepthPotential() { |
950 void updateDepthPotential() { |
837 depth[u_in] = depth[v_in] + 1; |
951 _depth[u_in] = _depth[v_in] + 1; |
838 potential[u_in] = forward[u_in] ? |
952 _potential[u_in] = _forward[u_in] ? |
839 potential[v_in] + cost[pred_edge[u_in]] : |
953 _potential[v_in] - _cost[_pred_edge[u_in]] : |
840 potential[v_in] - cost[pred_edge[u_in]]; |
954 _potential[v_in] + _cost[_pred_edge[u_in]]; |
841 |
955 |
842 Node u = thread[u_in], v; |
956 Node u = _thread[u_in], v; |
843 while (true) { |
957 while (true) { |
844 v = parent[u]; |
958 v = _parent[u]; |
845 if (v == INVALID) break; |
959 if (v == INVALID) break; |
846 depth[u] = depth[v] + 1; |
960 _depth[u] = _depth[v] + 1; |
847 potential[u] = forward[u] ? |
961 _potential[u] = _forward[u] ? |
848 potential[v] + cost[pred_edge[u]] : |
962 _potential[v] - _cost[_pred_edge[u]] : |
849 potential[v] - cost[pred_edge[u]]; |
963 _potential[v] + _cost[_pred_edge[u]]; |
850 if (depth[u] <= depth[v_in]) break; |
964 if (_depth[u] <= _depth[v_in]) break; |
851 u = thread[u]; |
965 u = _thread[u]; |
852 } |
966 } |
853 } |
967 } |
854 |
968 |
855 /// \brief Executes the algorithm. |
969 /// Executes the algorithm. |
|
970 bool start(PivotRuleEnum pivot_rule) { |
|
971 switch (pivot_rule) { |
|
972 case FIRST_ELIGIBLE_PIVOT: |
|
973 return start<FirstEligiblePivotRule>(); |
|
974 case BEST_ELIGIBLE_PIVOT: |
|
975 return start<BestEligiblePivotRule>(); |
|
976 case BLOCK_SEARCH_PIVOT: |
|
977 return start<BlockSearchPivotRule>(); |
|
978 case LIMITED_SEARCH_PIVOT: |
|
979 return start<LimitedSearchPivotRule>(); |
|
980 case CANDIDATE_LIST_PIVOT: |
|
981 return start<CandidateListPivotRule>(); |
|
982 case COMBINED_PIVOT: |
|
983 if ( countEdges(_graph) / countNodes(_graph) <= |
|
984 COMBINED_PIVOT_MAX_DEG ) |
|
985 return start<LimitedSearchPivotRule>(); |
|
986 else |
|
987 return start<BlockSearchPivotRule>(); |
|
988 } |
|
989 return false; |
|
990 } |
|
991 |
|
992 template<class PivotRuleImplementation> |
856 bool start() { |
993 bool start() { |
857 // Processing pivots |
994 PivotRuleImplementation pivot(*this); |
858 #ifdef _DEBUG_ITER_ |
995 |
859 int iter_num = 0; |
996 // Executing the network simplex algorithm |
860 #endif |
997 while (pivot.findEnteringEdge()) { |
861 #if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT) |
|
862 EdgeIt next_edge(graph); |
|
863 while (findEnteringEdge(next_edge)) |
|
864 #else |
|
865 while (findEnteringEdge()) |
|
866 #endif |
|
867 { |
|
868 join = findJoinNode(); |
998 join = findJoinNode(); |
869 bool change = findLeavingEdge(); |
999 bool change = findLeavingEdge(); |
870 changeFlows(change); |
1000 changeFlows(change); |
871 if (change) { |
1001 if (change) { |
872 updateThreadParent(); |
1002 updateThreadParent(); |
873 updatePredEdge(); |
1003 updatePredEdge(); |
874 updateDepthPotential(); |
1004 updateDepthPotential(); |
875 } |
1005 } |
876 #ifdef _DEBUG_ITER_ |
1006 } |
877 ++iter_num; |
1007 |
878 #endif |
1008 // Checking if the flow amount equals zero on all the artificial |
879 } |
1009 // edges |
880 |
1010 for (InEdgeIt e(_graph, _root); e != INVALID; ++e) |
881 #ifdef _DEBUG_ITER_ |
1011 if (_flow[e] > 0) return false; |
882 std::cout << "Network Simplex algorithm finished. " << iter_num |
1012 for (OutEdgeIt e(_graph, _root); e != INVALID; ++e) |
883 << " pivot iterations performed." << std::endl; |
1013 if (_flow[e] > 0) return false; |
884 #endif |
1014 |
885 |
1015 // Copying flow values to _flow_result |
886 // Checking if the flow amount equals zero on all the |
1016 if (_lower) { |
887 // artificial edges |
1017 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
888 for (InEdgeIt e(graph, root); e != INVALID; ++e) |
1018 _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]]; |
889 if (flow[e] > 0) return false; |
|
890 for (OutEdgeIt e(graph, root); e != INVALID; ++e) |
|
891 if (flow[e] > 0) return false; |
|
892 |
|
893 // Copying flow values to flow_result |
|
894 if (lower) { |
|
895 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) |
|
896 flow_result[e] = (*lower)[e] + flow[edge_ref[e]]; |
|
897 } else { |
1019 } else { |
898 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) |
1020 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
899 flow_result[e] = flow[edge_ref[e]]; |
1021 _flow_result[e] = _flow[_edge_ref[e]]; |
900 } |
1022 } |
901 // Copying potential values to potential_result |
1023 // Copying potential values to _potential_result |
902 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
1024 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
903 potential_result[n] = potential[node_ref[n]]; |
1025 _potential_result[n] = _potential[_node_ref[n]]; |
904 |
1026 |
905 return true; |
1027 return true; |
906 } |
1028 } |
907 |
1029 |
908 }; //class NetworkSimplex |
1030 }; //class NetworkSimplex |