|
1 /* -*- C++ -*- |
|
2 * |
|
3 * This file is a part of LEMON, a generic C++ optimization library |
|
4 * |
|
5 * Copyright (C) 2003-2008 |
|
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
7 * (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
8 * |
|
9 * Permission to use, modify and distribute this software is granted |
|
10 * provided that this copyright notice appears in all copies. For |
|
11 * precise terms see the accompanying LICENSE file. |
|
12 * |
|
13 * This software is provided "AS IS" with no warranty of any kind, |
|
14 * express or implied, and with no claim as to its suitability for any |
|
15 * purpose. |
|
16 * |
|
17 */ |
|
18 |
|
19 #ifndef LEMON_COST_SCALING_H |
|
20 #define LEMON_COST_SCALING_H |
|
21 |
|
22 /// \ingroup min_cost_flow_algs |
|
23 /// \file |
|
24 /// \brief Cost scaling algorithm for finding a minimum cost flow. |
|
25 |
|
26 #include <vector> |
|
27 #include <deque> |
|
28 #include <limits> |
|
29 |
|
30 #include <lemon/core.h> |
|
31 #include <lemon/maps.h> |
|
32 #include <lemon/math.h> |
|
33 #include <lemon/adaptors.h> |
|
34 #include <lemon/circulation.h> |
|
35 #include <lemon/bellman_ford.h> |
|
36 |
|
37 namespace lemon { |
|
38 |
|
39 /// \addtogroup min_cost_flow_algs |
|
40 /// @{ |
|
41 |
|
42 /// \brief Implementation of the cost scaling algorithm for finding a |
|
43 /// minimum cost flow. |
|
44 /// |
|
45 /// \ref CostScaling implements the cost scaling algorithm performing |
|
46 /// augment/push and relabel operations for finding a minimum cost |
|
47 /// flow. |
|
48 /// |
|
49 /// \tparam Digraph The digraph type the algorithm runs on. |
|
50 /// \tparam LowerMap The type of the lower bound map. |
|
51 /// \tparam CapacityMap The type of the capacity (upper bound) map. |
|
52 /// \tparam CostMap The type of the cost (length) map. |
|
53 /// \tparam SupplyMap The type of the supply map. |
|
54 /// |
|
55 /// \warning |
|
56 /// - Arc capacities and costs should be \e non-negative \e integers. |
|
57 /// - Supply values should be \e signed \e integers. |
|
58 /// - The value types of the maps should be convertible to each other. |
|
59 /// - \c CostMap::Value must be signed type. |
|
60 /// |
|
61 /// \note Arc costs are multiplied with the number of nodes during |
|
62 /// the algorithm so overflow problems may arise more easily than with |
|
63 /// other minimum cost flow algorithms. |
|
64 /// If it is available, <tt>long long int</tt> type is used instead of |
|
65 /// <tt>long int</tt> in the inside computations. |
|
66 /// |
|
67 /// \author Peter Kovacs |
|
68 template < typename Digraph, |
|
69 typename LowerMap = typename Digraph::template ArcMap<int>, |
|
70 typename CapacityMap = typename Digraph::template ArcMap<int>, |
|
71 typename CostMap = typename Digraph::template ArcMap<int>, |
|
72 typename SupplyMap = typename Digraph::template NodeMap<int> > |
|
73 class CostScaling |
|
74 { |
|
75 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
76 |
|
77 typedef typename CapacityMap::Value Capacity; |
|
78 typedef typename CostMap::Value Cost; |
|
79 typedef typename SupplyMap::Value Supply; |
|
80 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap; |
|
81 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap; |
|
82 |
|
83 typedef ResidualDigraph< const Digraph, |
|
84 CapacityArcMap, CapacityArcMap > ResDigraph; |
|
85 typedef typename ResDigraph::Arc ResArc; |
|
86 |
|
87 #if defined __GNUC__ && !defined __STRICT_ANSI__ |
|
88 typedef long long int LCost; |
|
89 #else |
|
90 typedef long int LCost; |
|
91 #endif |
|
92 typedef typename Digraph::template ArcMap<LCost> LargeCostMap; |
|
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<LCost> PotentialMap; |
|
100 |
|
101 private: |
|
102 |
|
103 /// \brief Map adaptor class for handling residual arc costs. |
|
104 /// |
|
105 /// Map adaptor class for handling residual arc costs. |
|
106 template <typename Map> |
|
107 class ResidualCostMap : public MapBase<ResArc, typename Map::Value> |
|
108 { |
|
109 private: |
|
110 |
|
111 const Map &_cost_map; |
|
112 |
|
113 public: |
|
114 |
|
115 ///\e |
|
116 ResidualCostMap(const Map &cost_map) : |
|
117 _cost_map(cost_map) {} |
|
118 |
|
119 ///\e |
|
120 inline typename Map::Value operator[](const ResArc &e) const { |
|
121 return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e]; |
|
122 } |
|
123 |
|
124 }; //class ResidualCostMap |
|
125 |
|
126 /// \brief Map adaptor class for handling reduced arc costs. |
|
127 /// |
|
128 /// Map adaptor class for handling reduced arc costs. |
|
129 class ReducedCostMap : public MapBase<Arc, LCost> |
|
130 { |
|
131 private: |
|
132 |
|
133 const Digraph &_gr; |
|
134 const LargeCostMap &_cost_map; |
|
135 const PotentialMap &_pot_map; |
|
136 |
|
137 public: |
|
138 |
|
139 ///\e |
|
140 ReducedCostMap( const Digraph &gr, |
|
141 const LargeCostMap &cost_map, |
|
142 const PotentialMap &pot_map ) : |
|
143 _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {} |
|
144 |
|
145 ///\e |
|
146 inline LCost operator[](const Arc &e) const { |
|
147 return _cost_map[e] + _pot_map[_gr.source(e)] |
|
148 - _pot_map[_gr.target(e)]; |
|
149 } |
|
150 |
|
151 }; //class ReducedCostMap |
|
152 |
|
153 private: |
|
154 |
|
155 // The digraph the algorithm runs on |
|
156 const Digraph &_graph; |
|
157 // The original lower bound map |
|
158 const LowerMap *_lower; |
|
159 // The modified capacity map |
|
160 CapacityArcMap _capacity; |
|
161 // The original cost map |
|
162 const CostMap &_orig_cost; |
|
163 // The scaled cost map |
|
164 LargeCostMap _cost; |
|
165 // The modified supply map |
|
166 SupplyNodeMap _supply; |
|
167 bool _valid_supply; |
|
168 |
|
169 // Arc map of the current flow |
|
170 FlowMap *_flow; |
|
171 bool _local_flow; |
|
172 // Node map of the current potentials |
|
173 PotentialMap *_potential; |
|
174 bool _local_potential; |
|
175 |
|
176 // The residual cost map |
|
177 ResidualCostMap<LargeCostMap> _res_cost; |
|
178 // The residual digraph |
|
179 ResDigraph *_res_graph; |
|
180 // The reduced cost map |
|
181 ReducedCostMap *_red_cost; |
|
182 // The excess map |
|
183 SupplyNodeMap _excess; |
|
184 // The epsilon parameter used for cost scaling |
|
185 LCost _epsilon; |
|
186 // The scaling factor |
|
187 int _alpha; |
|
188 |
|
189 public: |
|
190 |
|
191 /// \brief General constructor (with lower bounds). |
|
192 /// |
|
193 /// General constructor (with lower bounds). |
|
194 /// |
|
195 /// \param digraph The digraph the algorithm runs on. |
|
196 /// \param lower The lower bounds of the arcs. |
|
197 /// \param capacity The capacities (upper bounds) of the arcs. |
|
198 /// \param cost The cost (length) values of the arcs. |
|
199 /// \param supply The supply values of the nodes (signed). |
|
200 CostScaling( const Digraph &digraph, |
|
201 const LowerMap &lower, |
|
202 const CapacityMap &capacity, |
|
203 const CostMap &cost, |
|
204 const SupplyMap &supply ) : |
|
205 _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost), |
|
206 _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false), |
|
207 _potential(NULL), _local_potential(false), _res_cost(_cost), |
|
208 _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) |
|
209 { |
|
210 // Check the sum of supply values |
|
211 Supply sum = 0; |
|
212 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; |
|
213 _valid_supply = sum == 0; |
|
214 |
|
215 for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e]; |
|
216 for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n]; |
|
217 |
|
218 // Remove non-zero lower bounds |
|
219 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
220 if (lower[e] != 0) { |
|
221 _capacity[e] -= lower[e]; |
|
222 _supply[_graph.source(e)] -= lower[e]; |
|
223 _supply[_graph.target(e)] += lower[e]; |
|
224 } |
|
225 } |
|
226 } |
|
227 /* |
|
228 /// \brief General constructor (without lower bounds). |
|
229 /// |
|
230 /// General constructor (without lower bounds). |
|
231 /// |
|
232 /// \param digraph The digraph the algorithm runs on. |
|
233 /// \param capacity The capacities (upper bounds) of the arcs. |
|
234 /// \param cost The cost (length) values of the arcs. |
|
235 /// \param supply The supply values of the nodes (signed). |
|
236 CostScaling( const Digraph &digraph, |
|
237 const CapacityMap &capacity, |
|
238 const CostMap &cost, |
|
239 const SupplyMap &supply ) : |
|
240 _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost), |
|
241 _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false), |
|
242 _potential(NULL), _local_potential(false), _res_cost(_cost), |
|
243 _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) |
|
244 { |
|
245 // Check the sum of supply values |
|
246 Supply sum = 0; |
|
247 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; |
|
248 _valid_supply = sum == 0; |
|
249 } |
|
250 |
|
251 /// \brief Simple constructor (with lower bounds). |
|
252 /// |
|
253 /// Simple constructor (with lower bounds). |
|
254 /// |
|
255 /// \param digraph The digraph the algorithm runs on. |
|
256 /// \param lower The lower bounds of the arcs. |
|
257 /// \param capacity The capacities (upper bounds) of the arcs. |
|
258 /// \param cost The cost (length) values of the arcs. |
|
259 /// \param s The source node. |
|
260 /// \param t The target node. |
|
261 /// \param flow_value The required amount of flow from node \c s |
|
262 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
|
263 CostScaling( const Digraph &digraph, |
|
264 const LowerMap &lower, |
|
265 const CapacityMap &capacity, |
|
266 const CostMap &cost, |
|
267 Node s, Node t, |
|
268 Supply flow_value ) : |
|
269 _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost), |
|
270 _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false), |
|
271 _potential(NULL), _local_potential(false), _res_cost(_cost), |
|
272 _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) |
|
273 { |
|
274 // Remove non-zero lower bounds |
|
275 _supply[s] = flow_value; |
|
276 _supply[t] = -flow_value; |
|
277 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
278 if (lower[e] != 0) { |
|
279 _capacity[e] -= lower[e]; |
|
280 _supply[_graph.source(e)] -= lower[e]; |
|
281 _supply[_graph.target(e)] += lower[e]; |
|
282 } |
|
283 } |
|
284 _valid_supply = true; |
|
285 } |
|
286 |
|
287 /// \brief Simple constructor (without lower bounds). |
|
288 /// |
|
289 /// Simple constructor (without lower bounds). |
|
290 /// |
|
291 /// \param digraph The digraph the algorithm runs on. |
|
292 /// \param capacity The capacities (upper bounds) of the arcs. |
|
293 /// \param cost The cost (length) values of the arcs. |
|
294 /// \param s The source node. |
|
295 /// \param t The target node. |
|
296 /// \param flow_value The required amount of flow from node \c s |
|
297 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
|
298 CostScaling( const Digraph &digraph, |
|
299 const CapacityMap &capacity, |
|
300 const CostMap &cost, |
|
301 Node s, Node t, |
|
302 Supply flow_value ) : |
|
303 _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost), |
|
304 _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false), |
|
305 _potential(NULL), _local_potential(false), _res_cost(_cost), |
|
306 _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) |
|
307 { |
|
308 _supply[s] = flow_value; |
|
309 _supply[t] = -flow_value; |
|
310 _valid_supply = true; |
|
311 } |
|
312 */ |
|
313 /// Destructor. |
|
314 ~CostScaling() { |
|
315 if (_local_flow) delete _flow; |
|
316 if (_local_potential) delete _potential; |
|
317 delete _res_graph; |
|
318 delete _red_cost; |
|
319 } |
|
320 |
|
321 /// \brief Set the flow map. |
|
322 /// |
|
323 /// Set the flow map. |
|
324 /// |
|
325 /// \return \c (*this) |
|
326 CostScaling& flowMap(FlowMap &map) { |
|
327 if (_local_flow) { |
|
328 delete _flow; |
|
329 _local_flow = false; |
|
330 } |
|
331 _flow = ↦ |
|
332 return *this; |
|
333 } |
|
334 |
|
335 /// \brief Set the potential map. |
|
336 /// |
|
337 /// Set the potential map. |
|
338 /// |
|
339 /// \return \c (*this) |
|
340 CostScaling& potentialMap(PotentialMap &map) { |
|
341 if (_local_potential) { |
|
342 delete _potential; |
|
343 _local_potential = false; |
|
344 } |
|
345 _potential = ↦ |
|
346 return *this; |
|
347 } |
|
348 |
|
349 /// \name Execution control |
|
350 |
|
351 /// @{ |
|
352 |
|
353 /// \brief Run the algorithm. |
|
354 /// |
|
355 /// Run the algorithm. |
|
356 /// |
|
357 /// \param partial_augment By default the algorithm performs |
|
358 /// partial augment and relabel operations in the cost scaling |
|
359 /// phases. Set this parameter to \c false for using local push and |
|
360 /// relabel operations instead. |
|
361 /// |
|
362 /// \return \c true if a feasible flow can be found. |
|
363 bool run(bool partial_augment = true) { |
|
364 if (partial_augment) { |
|
365 return init() && startPartialAugment(); |
|
366 } else { |
|
367 return init() && startPushRelabel(); |
|
368 } |
|
369 } |
|
370 |
|
371 /// @} |
|
372 |
|
373 /// \name Query Functions |
|
374 /// The result of the algorithm can be obtained using these |
|
375 /// functions.\n |
|
376 /// \ref lemon::CostScaling::run() "run()" must be called before |
|
377 /// using them. |
|
378 |
|
379 /// @{ |
|
380 |
|
381 /// \brief Return a const reference to the arc map storing the |
|
382 /// found flow. |
|
383 /// |
|
384 /// Return a const reference to the arc map storing the found flow. |
|
385 /// |
|
386 /// \pre \ref run() must be called before using this function. |
|
387 const FlowMap& flowMap() const { |
|
388 return *_flow; |
|
389 } |
|
390 |
|
391 /// \brief Return a const reference to the node map storing the |
|
392 /// found potentials (the dual solution). |
|
393 /// |
|
394 /// Return a const reference to the node map storing the found |
|
395 /// potentials (the dual solution). |
|
396 /// |
|
397 /// \pre \ref run() must be called before using this function. |
|
398 const PotentialMap& potentialMap() const { |
|
399 return *_potential; |
|
400 } |
|
401 |
|
402 /// \brief Return the flow on the given arc. |
|
403 /// |
|
404 /// Return the flow on the given arc. |
|
405 /// |
|
406 /// \pre \ref run() must be called before using this function. |
|
407 Capacity flow(const Arc& arc) const { |
|
408 return (*_flow)[arc]; |
|
409 } |
|
410 |
|
411 /// \brief Return the potential of the given node. |
|
412 /// |
|
413 /// Return the potential of the given node. |
|
414 /// |
|
415 /// \pre \ref run() must be called before using this function. |
|
416 Cost potential(const Node& node) const { |
|
417 return (*_potential)[node]; |
|
418 } |
|
419 |
|
420 /// \brief Return the total cost of the found flow. |
|
421 /// |
|
422 /// Return the total cost of the found flow. The complexity of the |
|
423 /// function is \f$ O(e) \f$. |
|
424 /// |
|
425 /// \pre \ref run() must be called before using this function. |
|
426 Cost totalCost() const { |
|
427 Cost c = 0; |
|
428 for (ArcIt e(_graph); e != INVALID; ++e) |
|
429 c += (*_flow)[e] * _orig_cost[e]; |
|
430 return c; |
|
431 } |
|
432 |
|
433 /// @} |
|
434 |
|
435 private: |
|
436 |
|
437 /// Initialize the algorithm. |
|
438 bool init() { |
|
439 if (!_valid_supply) return false; |
|
440 // The scaling factor |
|
441 _alpha = 8; |
|
442 |
|
443 // Initialize flow and potential maps |
|
444 if (!_flow) { |
|
445 _flow = new FlowMap(_graph); |
|
446 _local_flow = true; |
|
447 } |
|
448 if (!_potential) { |
|
449 _potential = new PotentialMap(_graph); |
|
450 _local_potential = true; |
|
451 } |
|
452 |
|
453 _red_cost = new ReducedCostMap(_graph, _cost, *_potential); |
|
454 _res_graph = new ResDigraph(_graph, _capacity, *_flow); |
|
455 |
|
456 // Initialize the scaled cost map and the epsilon parameter |
|
457 Cost max_cost = 0; |
|
458 int node_num = countNodes(_graph); |
|
459 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
460 _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha; |
|
461 if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e]; |
|
462 } |
|
463 _epsilon = max_cost * node_num; |
|
464 |
|
465 // Find a feasible flow using Circulation |
|
466 Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap, |
|
467 SupplyMap > |
|
468 circulation( _graph, constMap<Arc>(Capacity(0)), _capacity, |
|
469 _supply ); |
|
470 return circulation.flowMap(*_flow).run(); |
|
471 } |
|
472 |
|
473 /// Execute the algorithm performing partial augmentation and |
|
474 /// relabel operations. |
|
475 bool startPartialAugment() { |
|
476 // Paramters for heuristics |
|
477 // const int BF_HEURISTIC_EPSILON_BOUND = 1000; |
|
478 // const int BF_HEURISTIC_BOUND_FACTOR = 3; |
|
479 // Maximum augment path length |
|
480 const int MAX_PATH_LENGTH = 4; |
|
481 |
|
482 // Variables |
|
483 typename Digraph::template NodeMap<Arc> pred_arc(_graph); |
|
484 typename Digraph::template NodeMap<bool> forward(_graph); |
|
485 typename Digraph::template NodeMap<OutArcIt> next_out(_graph); |
|
486 typename Digraph::template NodeMap<InArcIt> next_in(_graph); |
|
487 typename Digraph::template NodeMap<bool> next_dir(_graph); |
|
488 std::deque<Node> active_nodes; |
|
489 std::vector<Node> path_nodes; |
|
490 |
|
491 // int node_num = countNodes(_graph); |
|
492 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
|
493 1 : _epsilon / _alpha ) |
|
494 { |
|
495 /* |
|
496 // "Early Termination" heuristic: use Bellman-Ford algorithm |
|
497 // to check if the current flow is optimal |
|
498 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { |
|
499 typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap; |
|
500 ShiftCostMap shift_cost(_res_cost, 1); |
|
501 BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost); |
|
502 bf.init(0); |
|
503 bool done = false; |
|
504 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); |
|
505 for (int i = 0; i < K && !done; ++i) |
|
506 done = bf.processNextWeakRound(); |
|
507 if (done) break; |
|
508 } |
|
509 */ |
|
510 // Saturate arcs not satisfying the optimality condition |
|
511 Capacity delta; |
|
512 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
513 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { |
|
514 delta = _capacity[e] - (*_flow)[e]; |
|
515 _excess[_graph.source(e)] -= delta; |
|
516 _excess[_graph.target(e)] += delta; |
|
517 (*_flow)[e] = _capacity[e]; |
|
518 } |
|
519 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { |
|
520 _excess[_graph.target(e)] -= (*_flow)[e]; |
|
521 _excess[_graph.source(e)] += (*_flow)[e]; |
|
522 (*_flow)[e] = 0; |
|
523 } |
|
524 } |
|
525 |
|
526 // Find active nodes (i.e. nodes with positive excess) |
|
527 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
528 if (_excess[n] > 0) active_nodes.push_back(n); |
|
529 } |
|
530 |
|
531 // Initialize the next arc maps |
|
532 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
533 next_out[n] = OutArcIt(_graph, n); |
|
534 next_in[n] = InArcIt(_graph, n); |
|
535 next_dir[n] = true; |
|
536 } |
|
537 |
|
538 // Perform partial augment and relabel operations |
|
539 while (active_nodes.size() > 0) { |
|
540 // Select an active node (FIFO selection) |
|
541 if (_excess[active_nodes[0]] <= 0) { |
|
542 active_nodes.pop_front(); |
|
543 continue; |
|
544 } |
|
545 Node start = active_nodes[0]; |
|
546 path_nodes.clear(); |
|
547 path_nodes.push_back(start); |
|
548 |
|
549 // Find an augmenting path from the start node |
|
550 Node u, tip = start; |
|
551 LCost min_red_cost; |
|
552 while ( _excess[tip] >= 0 && |
|
553 int(path_nodes.size()) <= MAX_PATH_LENGTH ) |
|
554 { |
|
555 if (next_dir[tip]) { |
|
556 for (OutArcIt e = next_out[tip]; e != INVALID; ++e) { |
|
557 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { |
|
558 u = _graph.target(e); |
|
559 pred_arc[u] = e; |
|
560 forward[u] = true; |
|
561 next_out[tip] = e; |
|
562 tip = u; |
|
563 path_nodes.push_back(tip); |
|
564 goto next_step; |
|
565 } |
|
566 } |
|
567 next_dir[tip] = false; |
|
568 } |
|
569 for (InArcIt e = next_in[tip]; e != INVALID; ++e) { |
|
570 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { |
|
571 u = _graph.source(e); |
|
572 pred_arc[u] = e; |
|
573 forward[u] = false; |
|
574 next_in[tip] = e; |
|
575 tip = u; |
|
576 path_nodes.push_back(tip); |
|
577 goto next_step; |
|
578 } |
|
579 } |
|
580 |
|
581 // Relabel tip node |
|
582 min_red_cost = std::numeric_limits<LCost>::max() / 2; |
|
583 for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) { |
|
584 if ( _capacity[oe] - (*_flow)[oe] > 0 && |
|
585 (*_red_cost)[oe] < min_red_cost ) |
|
586 min_red_cost = (*_red_cost)[oe]; |
|
587 } |
|
588 for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) { |
|
589 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost) |
|
590 min_red_cost = -(*_red_cost)[ie]; |
|
591 } |
|
592 (*_potential)[tip] -= min_red_cost + _epsilon; |
|
593 |
|
594 // Reset the next arc maps |
|
595 next_out[tip] = OutArcIt(_graph, tip); |
|
596 next_in[tip] = InArcIt(_graph, tip); |
|
597 next_dir[tip] = true; |
|
598 |
|
599 // Step back |
|
600 if (tip != start) { |
|
601 path_nodes.pop_back(); |
|
602 tip = path_nodes[path_nodes.size()-1]; |
|
603 } |
|
604 |
|
605 next_step: |
|
606 continue; |
|
607 } |
|
608 |
|
609 // Augment along the found path (as much flow as possible) |
|
610 Capacity delta; |
|
611 for (int i = 1; i < int(path_nodes.size()); ++i) { |
|
612 u = path_nodes[i]; |
|
613 delta = forward[u] ? |
|
614 _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] : |
|
615 (*_flow)[pred_arc[u]]; |
|
616 delta = std::min(delta, _excess[path_nodes[i-1]]); |
|
617 (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta; |
|
618 _excess[path_nodes[i-1]] -= delta; |
|
619 _excess[u] += delta; |
|
620 if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u); |
|
621 } |
|
622 } |
|
623 } |
|
624 |
|
625 // Compute node potentials for the original costs |
|
626 ResidualCostMap<CostMap> res_cost(_orig_cost); |
|
627 BellmanFord< ResDigraph, ResidualCostMap<CostMap> > |
|
628 bf(*_res_graph, res_cost); |
|
629 bf.init(0); bf.start(); |
|
630 for (NodeIt n(_graph); n != INVALID; ++n) |
|
631 (*_potential)[n] = bf.dist(n); |
|
632 |
|
633 // Handle non-zero lower bounds |
|
634 if (_lower) { |
|
635 for (ArcIt e(_graph); e != INVALID; ++e) |
|
636 (*_flow)[e] += (*_lower)[e]; |
|
637 } |
|
638 return true; |
|
639 } |
|
640 |
|
641 /// Execute the algorithm performing push and relabel operations. |
|
642 bool startPushRelabel() { |
|
643 // Paramters for heuristics |
|
644 // const int BF_HEURISTIC_EPSILON_BOUND = 1000; |
|
645 // const int BF_HEURISTIC_BOUND_FACTOR = 3; |
|
646 |
|
647 typename Digraph::template NodeMap<bool> hyper(_graph, false); |
|
648 typename Digraph::template NodeMap<Arc> pred_arc(_graph); |
|
649 typename Digraph::template NodeMap<bool> forward(_graph); |
|
650 typename Digraph::template NodeMap<OutArcIt> next_out(_graph); |
|
651 typename Digraph::template NodeMap<InArcIt> next_in(_graph); |
|
652 typename Digraph::template NodeMap<bool> next_dir(_graph); |
|
653 std::deque<Node> active_nodes; |
|
654 |
|
655 // int node_num = countNodes(_graph); |
|
656 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
|
657 1 : _epsilon / _alpha ) |
|
658 { |
|
659 /* |
|
660 // "Early Termination" heuristic: use Bellman-Ford algorithm |
|
661 // to check if the current flow is optimal |
|
662 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { |
|
663 typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap; |
|
664 ShiftCostMap shift_cost(_res_cost, 1); |
|
665 BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost); |
|
666 bf.init(0); |
|
667 bool done = false; |
|
668 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); |
|
669 for (int i = 0; i < K && !done; ++i) |
|
670 done = bf.processNextWeakRound(); |
|
671 if (done) break; |
|
672 } |
|
673 */ |
|
674 |
|
675 // Saturate arcs not satisfying the optimality condition |
|
676 Capacity delta; |
|
677 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
678 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { |
|
679 delta = _capacity[e] - (*_flow)[e]; |
|
680 _excess[_graph.source(e)] -= delta; |
|
681 _excess[_graph.target(e)] += delta; |
|
682 (*_flow)[e] = _capacity[e]; |
|
683 } |
|
684 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { |
|
685 _excess[_graph.target(e)] -= (*_flow)[e]; |
|
686 _excess[_graph.source(e)] += (*_flow)[e]; |
|
687 (*_flow)[e] = 0; |
|
688 } |
|
689 } |
|
690 |
|
691 // Find active nodes (i.e. nodes with positive excess) |
|
692 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
693 if (_excess[n] > 0) active_nodes.push_back(n); |
|
694 } |
|
695 |
|
696 // Initialize the next arc maps |
|
697 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
698 next_out[n] = OutArcIt(_graph, n); |
|
699 next_in[n] = InArcIt(_graph, n); |
|
700 next_dir[n] = true; |
|
701 } |
|
702 |
|
703 // Perform push and relabel operations |
|
704 while (active_nodes.size() > 0) { |
|
705 // Select an active node (FIFO selection) |
|
706 Node n = active_nodes[0], t; |
|
707 bool relabel_enabled = true; |
|
708 |
|
709 // Perform push operations if there are admissible arcs |
|
710 if (_excess[n] > 0 && next_dir[n]) { |
|
711 OutArcIt e = next_out[n]; |
|
712 for ( ; e != INVALID; ++e) { |
|
713 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { |
|
714 delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]); |
|
715 t = _graph.target(e); |
|
716 |
|
717 // Push-look-ahead heuristic |
|
718 Capacity ahead = -_excess[t]; |
|
719 for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) { |
|
720 if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0) |
|
721 ahead += _capacity[oe] - (*_flow)[oe]; |
|
722 } |
|
723 for (InArcIt ie(_graph, t); ie != INVALID; ++ie) { |
|
724 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0) |
|
725 ahead += (*_flow)[ie]; |
|
726 } |
|
727 if (ahead < 0) ahead = 0; |
|
728 |
|
729 // Push flow along the arc |
|
730 if (ahead < delta) { |
|
731 (*_flow)[e] += ahead; |
|
732 _excess[n] -= ahead; |
|
733 _excess[t] += ahead; |
|
734 active_nodes.push_front(t); |
|
735 hyper[t] = true; |
|
736 relabel_enabled = false; |
|
737 break; |
|
738 } else { |
|
739 (*_flow)[e] += delta; |
|
740 _excess[n] -= delta; |
|
741 _excess[t] += delta; |
|
742 if (_excess[t] > 0 && _excess[t] <= delta) |
|
743 active_nodes.push_back(t); |
|
744 } |
|
745 |
|
746 if (_excess[n] == 0) break; |
|
747 } |
|
748 } |
|
749 if (e != INVALID) { |
|
750 next_out[n] = e; |
|
751 } else { |
|
752 next_dir[n] = false; |
|
753 } |
|
754 } |
|
755 |
|
756 if (_excess[n] > 0 && !next_dir[n]) { |
|
757 InArcIt e = next_in[n]; |
|
758 for ( ; e != INVALID; ++e) { |
|
759 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { |
|
760 delta = std::min((*_flow)[e], _excess[n]); |
|
761 t = _graph.source(e); |
|
762 |
|
763 // Push-look-ahead heuristic |
|
764 Capacity ahead = -_excess[t]; |
|
765 for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) { |
|
766 if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0) |
|
767 ahead += _capacity[oe] - (*_flow)[oe]; |
|
768 } |
|
769 for (InArcIt ie(_graph, t); ie != INVALID; ++ie) { |
|
770 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0) |
|
771 ahead += (*_flow)[ie]; |
|
772 } |
|
773 if (ahead < 0) ahead = 0; |
|
774 |
|
775 // Push flow along the arc |
|
776 if (ahead < delta) { |
|
777 (*_flow)[e] -= ahead; |
|
778 _excess[n] -= ahead; |
|
779 _excess[t] += ahead; |
|
780 active_nodes.push_front(t); |
|
781 hyper[t] = true; |
|
782 relabel_enabled = false; |
|
783 break; |
|
784 } else { |
|
785 (*_flow)[e] -= delta; |
|
786 _excess[n] -= delta; |
|
787 _excess[t] += delta; |
|
788 if (_excess[t] > 0 && _excess[t] <= delta) |
|
789 active_nodes.push_back(t); |
|
790 } |
|
791 |
|
792 if (_excess[n] == 0) break; |
|
793 } |
|
794 } |
|
795 next_in[n] = e; |
|
796 } |
|
797 |
|
798 // Relabel the node if it is still active (or hyper) |
|
799 if (relabel_enabled && (_excess[n] > 0 || hyper[n])) { |
|
800 LCost min_red_cost = std::numeric_limits<LCost>::max() / 2; |
|
801 for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) { |
|
802 if ( _capacity[oe] - (*_flow)[oe] > 0 && |
|
803 (*_red_cost)[oe] < min_red_cost ) |
|
804 min_red_cost = (*_red_cost)[oe]; |
|
805 } |
|
806 for (InArcIt ie(_graph, n); ie != INVALID; ++ie) { |
|
807 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost) |
|
808 min_red_cost = -(*_red_cost)[ie]; |
|
809 } |
|
810 (*_potential)[n] -= min_red_cost + _epsilon; |
|
811 hyper[n] = false; |
|
812 |
|
813 // Reset the next arc maps |
|
814 next_out[n] = OutArcIt(_graph, n); |
|
815 next_in[n] = InArcIt(_graph, n); |
|
816 next_dir[n] = true; |
|
817 } |
|
818 |
|
819 // Remove nodes that are not active nor hyper |
|
820 while ( active_nodes.size() > 0 && |
|
821 _excess[active_nodes[0]] <= 0 && |
|
822 !hyper[active_nodes[0]] ) { |
|
823 active_nodes.pop_front(); |
|
824 } |
|
825 } |
|
826 } |
|
827 |
|
828 // Compute node potentials for the original costs |
|
829 ResidualCostMap<CostMap> res_cost(_orig_cost); |
|
830 BellmanFord< ResDigraph, ResidualCostMap<CostMap> > |
|
831 bf(*_res_graph, res_cost); |
|
832 bf.init(0); bf.start(); |
|
833 for (NodeIt n(_graph); n != INVALID; ++n) |
|
834 (*_potential)[n] = bf.dist(n); |
|
835 |
|
836 // Handle non-zero lower bounds |
|
837 if (_lower) { |
|
838 for (ArcIt e(_graph); e != INVALID; ++e) |
|
839 (*_flow)[e] += (*_lower)[e]; |
|
840 } |
|
841 return true; |
|
842 } |
|
843 |
|
844 }; //class CostScaling |
|
845 |
|
846 ///@} |
|
847 |
|
848 } //namespace lemon |
|
849 |
|
850 #endif //LEMON_COST_SCALING_H |