17 */ |
17 */ |
18 |
18 |
19 #ifndef LEMON_CAPACITY_SCALING_H |
19 #ifndef LEMON_CAPACITY_SCALING_H |
20 #define LEMON_CAPACITY_SCALING_H |
20 #define LEMON_CAPACITY_SCALING_H |
21 |
21 |
22 /// \ingroup min_cost_flow |
22 /// \ingroup min_cost_flow_algs |
23 /// |
23 /// |
24 /// \file |
24 /// \file |
25 /// \brief Capacity scaling algorithm for finding a minimum cost flow. |
25 /// \brief Capacity Scaling algorithm for finding a minimum cost flow. |
26 |
26 |
27 #include <vector> |
27 #include <vector> |
|
28 #include <limits> |
|
29 #include <lemon/core.h> |
28 #include <lemon/bin_heap.h> |
30 #include <lemon/bin_heap.h> |
29 |
31 |
30 namespace lemon { |
32 namespace lemon { |
31 |
33 |
32 /// \addtogroup min_cost_flow |
34 /// \addtogroup min_cost_flow_algs |
33 /// @{ |
35 /// @{ |
34 |
36 |
35 /// \brief Implementation of the capacity scaling algorithm for |
37 /// \brief Implementation of the Capacity Scaling algorithm for |
36 /// finding a minimum cost flow. |
38 /// finding a \ref min_cost_flow "minimum cost flow". |
37 /// |
39 /// |
38 /// \ref CapacityScaling implements the capacity scaling version |
40 /// \ref CapacityScaling implements the capacity scaling version |
39 /// of the successive shortest path algorithm for finding a minimum |
41 /// of the successive shortest path algorithm for finding a |
40 /// cost flow. |
42 /// \ref min_cost_flow "minimum cost flow". It is an efficient dual |
|
43 /// solution method. |
41 /// |
44 /// |
42 /// \tparam Digraph The digraph type the algorithm runs on. |
45 /// Most of the parameters of the problem (except for the digraph) |
43 /// \tparam LowerMap The type of the lower bound map. |
46 /// can be given using separate functions, and the algorithm can be |
44 /// \tparam CapacityMap The type of the capacity (upper bound) map. |
47 /// executed using the \ref run() function. If some parameters are not |
45 /// \tparam CostMap The type of the cost (length) map. |
48 /// specified, then default values will be used. |
46 /// \tparam SupplyMap The type of the supply map. |
|
47 /// |
49 /// |
48 /// \warning |
50 /// \tparam GR The digraph type the algorithm runs on. |
49 /// - Arc capacities and costs should be \e non-negative \e integers. |
51 /// \tparam V The value type used for flow amounts, capacity bounds |
50 /// - Supply values should be \e signed \e integers. |
52 /// and supply values in the algorithm. By default it is \c int. |
51 /// - The value types of the maps should be convertible to each other. |
53 /// \tparam C The value type used for costs and potentials in the |
52 /// - \c CostMap::Value must be signed type. |
54 /// algorithm. By default it is the same as \c V. |
53 /// |
55 /// |
54 /// \author Peter Kovacs |
56 /// \warning Both value types must be signed and all input data must |
55 template < typename Digraph, |
57 /// be integer. |
56 typename LowerMap = typename Digraph::template ArcMap<int>, |
58 /// \warning This algorithm does not support negative costs for such |
57 typename CapacityMap = typename Digraph::template ArcMap<int>, |
59 /// arcs that have infinite upper bound. |
58 typename CostMap = typename Digraph::template ArcMap<int>, |
60 template <typename GR, typename V = int, typename C = V> |
59 typename SupplyMap = typename Digraph::template NodeMap<int> > |
|
60 class CapacityScaling |
61 class CapacityScaling |
61 { |
62 { |
62 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
63 |
|
64 typedef typename CapacityMap::Value Capacity; |
|
65 typedef typename CostMap::Value Cost; |
|
66 typedef typename SupplyMap::Value Supply; |
|
67 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap; |
|
68 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap; |
|
69 typedef typename Digraph::template NodeMap<Arc> PredMap; |
|
70 |
|
71 public: |
63 public: |
72 |
64 |
73 /// The type of the flow map. |
65 /// The type of the flow amounts, capacity bounds and supply values |
74 typedef typename Digraph::template ArcMap<Capacity> FlowMap; |
66 typedef V Value; |
75 /// The type of the potential map. |
67 /// The type of the arc costs |
76 typedef typename Digraph::template NodeMap<Cost> PotentialMap; |
68 typedef C Cost; |
77 |
69 |
|
70 public: |
|
71 |
|
72 /// \brief Problem type constants for the \c run() function. |
|
73 /// |
|
74 /// Enum type containing the problem type constants that can be |
|
75 /// returned by the \ref run() function of the algorithm. |
|
76 enum ProblemType { |
|
77 /// The problem has no feasible solution (flow). |
|
78 INFEASIBLE, |
|
79 /// The problem has optimal solution (i.e. it is feasible and |
|
80 /// bounded), and the algorithm has found optimal flow and node |
|
81 /// potentials (primal and dual solutions). |
|
82 OPTIMAL, |
|
83 /// The digraph contains an arc of negative cost and infinite |
|
84 /// upper bound. It means that the objective function is unbounded |
|
85 /// on that arc, however note that it could actually be bounded |
|
86 /// over the feasible flows, but this algroithm cannot handle |
|
87 /// these cases. |
|
88 UNBOUNDED |
|
89 }; |
|
90 |
78 private: |
91 private: |
79 |
92 |
80 /// \brief Special implementation of the \ref Dijkstra algorithm |
93 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
81 /// for finding shortest paths in the residual network. |
94 |
82 /// |
95 typedef std::vector<Arc> ArcVector; |
83 /// \ref ResidualDijkstra is a special implementation of the |
96 typedef std::vector<Node> NodeVector; |
84 /// \ref Dijkstra algorithm for finding shortest paths in the |
97 typedef std::vector<int> IntVector; |
85 /// residual network of the digraph with respect to the reduced arc |
98 typedef std::vector<bool> BoolVector; |
86 /// costs and modifying the node potentials according to the |
99 typedef std::vector<Value> ValueVector; |
87 /// distance of the nodes. |
100 typedef std::vector<Cost> CostVector; |
|
101 |
|
102 private: |
|
103 |
|
104 // Data related to the underlying digraph |
|
105 const GR &_graph; |
|
106 int _node_num; |
|
107 int _arc_num; |
|
108 int _res_arc_num; |
|
109 int _root; |
|
110 |
|
111 // Parameters of the problem |
|
112 bool _have_lower; |
|
113 Value _sum_supply; |
|
114 |
|
115 // Data structures for storing the digraph |
|
116 IntNodeMap _node_id; |
|
117 IntArcMap _arc_idf; |
|
118 IntArcMap _arc_idb; |
|
119 IntVector _first_out; |
|
120 BoolVector _forward; |
|
121 IntVector _source; |
|
122 IntVector _target; |
|
123 IntVector _reverse; |
|
124 |
|
125 // Node and arc data |
|
126 ValueVector _lower; |
|
127 ValueVector _upper; |
|
128 CostVector _cost; |
|
129 ValueVector _supply; |
|
130 |
|
131 ValueVector _res_cap; |
|
132 CostVector _pi; |
|
133 ValueVector _excess; |
|
134 IntVector _excess_nodes; |
|
135 IntVector _deficit_nodes; |
|
136 |
|
137 Value _delta; |
|
138 int _phase_num; |
|
139 IntVector _pred; |
|
140 |
|
141 public: |
|
142 |
|
143 /// \brief Constant for infinite upper bounds (capacities). |
|
144 /// |
|
145 /// Constant for infinite upper bounds (capacities). |
|
146 /// It is \c std::numeric_limits<Value>::infinity() if available, |
|
147 /// \c std::numeric_limits<Value>::max() otherwise. |
|
148 const Value INF; |
|
149 |
|
150 private: |
|
151 |
|
152 // Special implementation of the Dijkstra algorithm for finding |
|
153 // shortest paths in the residual network of the digraph with |
|
154 // respect to the reduced arc costs and modifying the node |
|
155 // potentials according to the found distance labels. |
88 class ResidualDijkstra |
156 class ResidualDijkstra |
89 { |
157 { |
90 typedef typename Digraph::template NodeMap<int> HeapCrossRef; |
158 typedef RangeMap<int> HeapCrossRef; |
91 typedef BinHeap<Cost, HeapCrossRef> Heap; |
159 typedef BinHeap<Cost, HeapCrossRef> Heap; |
92 |
160 |
93 private: |
161 private: |
94 |
162 |
95 // The digraph the algorithm runs on |
163 int _node_num; |
96 const Digraph &_graph; |
164 const IntVector &_first_out; |
97 |
165 const IntVector &_target; |
98 // The main maps |
166 const CostVector &_cost; |
99 const FlowMap &_flow; |
167 const ValueVector &_res_cap; |
100 const CapacityArcMap &_res_cap; |
168 const ValueVector &_excess; |
101 const CostMap &_cost; |
169 CostVector &_pi; |
102 const SupplyNodeMap &_excess; |
170 IntVector &_pred; |
103 PotentialMap &_potential; |
171 |
104 |
172 IntVector _proc_nodes; |
105 // The distance map |
173 CostVector _dist; |
106 PotentialMap _dist; |
174 |
107 // The pred arc map |
|
108 PredMap &_pred; |
|
109 // The processed (i.e. permanently labeled) nodes |
|
110 std::vector<Node> _proc_nodes; |
|
111 |
|
112 public: |
175 public: |
113 |
176 |
114 /// Constructor. |
177 ResidualDijkstra(CapacityScaling& cs) : |
115 ResidualDijkstra( const Digraph &digraph, |
178 _node_num(cs._node_num), _first_out(cs._first_out), |
116 const FlowMap &flow, |
179 _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap), |
117 const CapacityArcMap &res_cap, |
180 _excess(cs._excess), _pi(cs._pi), _pred(cs._pred), |
118 const CostMap &cost, |
181 _dist(cs._node_num) |
119 const SupplyMap &excess, |
|
120 PotentialMap &potential, |
|
121 PredMap &pred ) : |
|
122 _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost), |
|
123 _excess(excess), _potential(potential), _dist(digraph), |
|
124 _pred(pred) |
|
125 {} |
182 {} |
126 |
183 |
127 /// Run the algorithm from the given source node. |
184 int run(int s, Value delta = 1) { |
128 Node run(Node s, Capacity delta = 1) { |
185 HeapCrossRef heap_cross_ref(_node_num, Heap::PRE_HEAP); |
129 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); |
|
130 Heap heap(heap_cross_ref); |
186 Heap heap(heap_cross_ref); |
131 heap.push(s, 0); |
187 heap.push(s, 0); |
132 _pred[s] = INVALID; |
188 _pred[s] = -1; |
133 _proc_nodes.clear(); |
189 _proc_nodes.clear(); |
134 |
190 |
135 // Processing nodes |
191 // Process nodes |
136 while (!heap.empty() && _excess[heap.top()] > -delta) { |
192 while (!heap.empty() && _excess[heap.top()] > -delta) { |
137 Node u = heap.top(), v; |
193 int u = heap.top(), v; |
138 Cost d = heap.prio() + _potential[u], nd; |
194 Cost d = heap.prio() + _pi[u], dn; |
139 _dist[u] = heap.prio(); |
195 _dist[u] = heap.prio(); |
|
196 _proc_nodes.push_back(u); |
140 heap.pop(); |
197 heap.pop(); |
141 _proc_nodes.push_back(u); |
198 |
142 |
199 // Traverse outgoing residual arcs |
143 // Traversing outgoing arcs |
200 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { |
144 for (OutArcIt e(_graph, u); e != INVALID; ++e) { |
201 if (_res_cap[a] < delta) continue; |
145 if (_res_cap[e] >= delta) { |
202 v = _target[a]; |
146 v = _graph.target(e); |
203 switch (heap.state(v)) { |
147 switch(heap.state(v)) { |
|
148 case Heap::PRE_HEAP: |
204 case Heap::PRE_HEAP: |
149 heap.push(v, d + _cost[e] - _potential[v]); |
205 heap.push(v, d + _cost[a] - _pi[v]); |
150 _pred[v] = e; |
206 _pred[v] = a; |
151 break; |
207 break; |
152 case Heap::IN_HEAP: |
208 case Heap::IN_HEAP: |
153 nd = d + _cost[e] - _potential[v]; |
209 dn = d + _cost[a] - _pi[v]; |
154 if (nd < heap[v]) { |
210 if (dn < heap[v]) { |
155 heap.decrease(v, nd); |
211 heap.decrease(v, dn); |
156 _pred[v] = e; |
212 _pred[v] = a; |
157 } |
213 } |
158 break; |
214 break; |
159 case Heap::POST_HEAP: |
215 case Heap::POST_HEAP: |
160 break; |
216 break; |
161 } |
|
162 } |
217 } |
163 } |
218 } |
164 |
219 } |
165 // Traversing incoming arcs |
220 if (heap.empty()) return -1; |
166 for (InArcIt e(_graph, u); e != INVALID; ++e) { |
221 |
167 if (_flow[e] >= delta) { |
222 // Update potentials of processed nodes |
168 v = _graph.source(e); |
223 int t = heap.top(); |
169 switch(heap.state(v)) { |
224 Cost dt = heap.prio(); |
170 case Heap::PRE_HEAP: |
225 for (int i = 0; i < int(_proc_nodes.size()); ++i) { |
171 heap.push(v, d - _cost[e] - _potential[v]); |
226 _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt; |
172 _pred[v] = e; |
227 } |
173 break; |
|
174 case Heap::IN_HEAP: |
|
175 nd = d - _cost[e] - _potential[v]; |
|
176 if (nd < heap[v]) { |
|
177 heap.decrease(v, nd); |
|
178 _pred[v] = e; |
|
179 } |
|
180 break; |
|
181 case Heap::POST_HEAP: |
|
182 break; |
|
183 } |
|
184 } |
|
185 } |
|
186 } |
|
187 if (heap.empty()) return INVALID; |
|
188 |
|
189 // Updating potentials of processed nodes |
|
190 Node t = heap.top(); |
|
191 Cost t_dist = heap.prio(); |
|
192 for (int i = 0; i < int(_proc_nodes.size()); ++i) |
|
193 _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; |
|
194 |
228 |
195 return t; |
229 return t; |
196 } |
230 } |
197 |
231 |
198 }; //class ResidualDijkstra |
232 }; //class ResidualDijkstra |
199 |
233 |
200 private: |
|
201 |
|
202 // The digraph the algorithm runs on |
|
203 const Digraph &_graph; |
|
204 // The original lower bound map |
|
205 const LowerMap *_lower; |
|
206 // The modified capacity map |
|
207 CapacityArcMap _capacity; |
|
208 // The original cost map |
|
209 const CostMap &_cost; |
|
210 // The modified supply map |
|
211 SupplyNodeMap _supply; |
|
212 bool _valid_supply; |
|
213 |
|
214 // Arc map of the current flow |
|
215 FlowMap *_flow; |
|
216 bool _local_flow; |
|
217 // Node map of the current potentials |
|
218 PotentialMap *_potential; |
|
219 bool _local_potential; |
|
220 |
|
221 // The residual capacity map |
|
222 CapacityArcMap _res_cap; |
|
223 // The excess map |
|
224 SupplyNodeMap _excess; |
|
225 // The excess nodes (i.e. nodes with positive excess) |
|
226 std::vector<Node> _excess_nodes; |
|
227 // The deficit nodes (i.e. nodes with negative excess) |
|
228 std::vector<Node> _deficit_nodes; |
|
229 |
|
230 // The delta parameter used for capacity scaling |
|
231 Capacity _delta; |
|
232 // The maximum number of phases |
|
233 int _phase_num; |
|
234 |
|
235 // The pred arc map |
|
236 PredMap _pred; |
|
237 // Implementation of the Dijkstra algorithm for finding augmenting |
|
238 // shortest paths in the residual network |
|
239 ResidualDijkstra *_dijkstra; |
|
240 |
|
241 public: |
234 public: |
242 |
235 |
243 /// \brief General constructor (with lower bounds). |
236 /// \brief Constructor. |
244 /// |
237 /// |
245 /// General constructor (with lower bounds). |
238 /// The constructor of the class. |
246 /// |
239 /// |
247 /// \param digraph The digraph the algorithm runs on. |
240 /// \param graph The digraph the algorithm runs on. |
248 /// \param lower The lower bounds of the arcs. |
241 CapacityScaling(const GR& graph) : |
249 /// \param capacity The capacities (upper bounds) of the arcs. |
242 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
250 /// \param cost The cost (length) values of the arcs. |
243 INF(std::numeric_limits<Value>::has_infinity ? |
251 /// \param supply The supply values of the nodes (signed). |
244 std::numeric_limits<Value>::infinity() : |
252 CapacityScaling( const Digraph &digraph, |
245 std::numeric_limits<Value>::max()) |
253 const LowerMap &lower, |
|
254 const CapacityMap &capacity, |
|
255 const CostMap &cost, |
|
256 const SupplyMap &supply ) : |
|
257 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost), |
|
258 _supply(digraph), _flow(NULL), _local_flow(false), |
|
259 _potential(NULL), _local_potential(false), |
|
260 _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL) |
|
261 { |
246 { |
262 Supply sum = 0; |
247 // Check the value types |
|
248 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
|
249 "The flow type of CapacityScaling must be signed"); |
|
250 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
|
251 "The cost type of CapacityScaling must be signed"); |
|
252 |
|
253 // Resize vectors |
|
254 _node_num = countNodes(_graph); |
|
255 _arc_num = countArcs(_graph); |
|
256 _res_arc_num = 2 * (_arc_num + _node_num); |
|
257 _root = _node_num; |
|
258 ++_node_num; |
|
259 |
|
260 _first_out.resize(_node_num + 1); |
|
261 _forward.resize(_res_arc_num); |
|
262 _source.resize(_res_arc_num); |
|
263 _target.resize(_res_arc_num); |
|
264 _reverse.resize(_res_arc_num); |
|
265 |
|
266 _lower.resize(_res_arc_num); |
|
267 _upper.resize(_res_arc_num); |
|
268 _cost.resize(_res_arc_num); |
|
269 _supply.resize(_node_num); |
|
270 |
|
271 _res_cap.resize(_res_arc_num); |
|
272 _pi.resize(_node_num); |
|
273 _excess.resize(_node_num); |
|
274 _pred.resize(_node_num); |
|
275 |
|
276 // Copy the graph |
|
277 int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1; |
|
278 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
279 _node_id[n] = i; |
|
280 } |
|
281 i = 0; |
|
282 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
283 _first_out[i] = j; |
|
284 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
285 _arc_idf[a] = j; |
|
286 _forward[j] = true; |
|
287 _source[j] = i; |
|
288 _target[j] = _node_id[_graph.runningNode(a)]; |
|
289 } |
|
290 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
291 _arc_idb[a] = j; |
|
292 _forward[j] = false; |
|
293 _source[j] = i; |
|
294 _target[j] = _node_id[_graph.runningNode(a)]; |
|
295 } |
|
296 _forward[j] = false; |
|
297 _source[j] = i; |
|
298 _target[j] = _root; |
|
299 _reverse[j] = k; |
|
300 _forward[k] = true; |
|
301 _source[k] = _root; |
|
302 _target[k] = i; |
|
303 _reverse[k] = j; |
|
304 ++j; ++k; |
|
305 } |
|
306 _first_out[i] = j; |
|
307 _first_out[_node_num] = k; |
|
308 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
309 int fi = _arc_idf[a]; |
|
310 int bi = _arc_idb[a]; |
|
311 _reverse[fi] = bi; |
|
312 _reverse[bi] = fi; |
|
313 } |
|
314 |
|
315 // Reset parameters |
|
316 reset(); |
|
317 } |
|
318 |
|
319 /// \name Parameters |
|
320 /// The parameters of the algorithm can be specified using these |
|
321 /// functions. |
|
322 |
|
323 /// @{ |
|
324 |
|
325 /// \brief Set the lower bounds on the arcs. |
|
326 /// |
|
327 /// This function sets the lower bounds on the arcs. |
|
328 /// If it is not used before calling \ref run(), the lower bounds |
|
329 /// will be set to zero on all arcs. |
|
330 /// |
|
331 /// \param map An arc map storing the lower bounds. |
|
332 /// Its \c Value type must be convertible to the \c Value type |
|
333 /// of the algorithm. |
|
334 /// |
|
335 /// \return <tt>(*this)</tt> |
|
336 template <typename LowerMap> |
|
337 CapacityScaling& lowerMap(const LowerMap& map) { |
|
338 _have_lower = true; |
|
339 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
340 _lower[_arc_idf[a]] = map[a]; |
|
341 _lower[_arc_idb[a]] = map[a]; |
|
342 } |
|
343 return *this; |
|
344 } |
|
345 |
|
346 /// \brief Set the upper bounds (capacities) on the arcs. |
|
347 /// |
|
348 /// This function sets the upper bounds (capacities) on the arcs. |
|
349 /// If it is not used before calling \ref run(), the upper bounds |
|
350 /// will be set to \ref INF on all arcs (i.e. the flow value will be |
|
351 /// unbounded from above on each arc). |
|
352 /// |
|
353 /// \param map An arc map storing the upper bounds. |
|
354 /// Its \c Value type must be convertible to the \c Value type |
|
355 /// of the algorithm. |
|
356 /// |
|
357 /// \return <tt>(*this)</tt> |
|
358 template<typename UpperMap> |
|
359 CapacityScaling& upperMap(const UpperMap& map) { |
|
360 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
361 _upper[_arc_idf[a]] = map[a]; |
|
362 } |
|
363 return *this; |
|
364 } |
|
365 |
|
366 /// \brief Set the costs of the arcs. |
|
367 /// |
|
368 /// This function sets the costs of the arcs. |
|
369 /// If it is not used before calling \ref run(), the costs |
|
370 /// will be set to \c 1 on all arcs. |
|
371 /// |
|
372 /// \param map An arc map storing the costs. |
|
373 /// Its \c Value type must be convertible to the \c Cost type |
|
374 /// of the algorithm. |
|
375 /// |
|
376 /// \return <tt>(*this)</tt> |
|
377 template<typename CostMap> |
|
378 CapacityScaling& costMap(const CostMap& map) { |
|
379 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
380 _cost[_arc_idf[a]] = map[a]; |
|
381 _cost[_arc_idb[a]] = -map[a]; |
|
382 } |
|
383 return *this; |
|
384 } |
|
385 |
|
386 /// \brief Set the supply values of the nodes. |
|
387 /// |
|
388 /// This function sets the supply values of the nodes. |
|
389 /// If neither this function nor \ref stSupply() is used before |
|
390 /// calling \ref run(), the supply of each node will be set to zero. |
|
391 /// |
|
392 /// \param map A node map storing the supply values. |
|
393 /// Its \c Value type must be convertible to the \c Value type |
|
394 /// of the algorithm. |
|
395 /// |
|
396 /// \return <tt>(*this)</tt> |
|
397 template<typename SupplyMap> |
|
398 CapacityScaling& supplyMap(const SupplyMap& map) { |
263 for (NodeIt n(_graph); n != INVALID; ++n) { |
399 for (NodeIt n(_graph); n != INVALID; ++n) { |
264 _supply[n] = supply[n]; |
400 _supply[_node_id[n]] = map[n]; |
265 _excess[n] = supply[n]; |
401 } |
266 sum += supply[n]; |
402 return *this; |
267 } |
403 } |
268 _valid_supply = sum == 0; |
404 |
269 for (ArcIt a(_graph); a != INVALID; ++a) { |
405 /// \brief Set single source and target nodes and a supply value. |
270 _capacity[a] = capacity[a]; |
406 /// |
271 _res_cap[a] = capacity[a]; |
407 /// This function sets a single source node and a single target node |
272 } |
408 /// and the required flow value. |
273 |
409 /// If neither this function nor \ref supplyMap() is used before |
274 // Remove non-zero lower bounds |
410 /// calling \ref run(), the supply of each node will be set to zero. |
275 typename LowerMap::Value lcap; |
411 /// |
276 for (ArcIt e(_graph); e != INVALID; ++e) { |
412 /// Using this function has the same effect as using \ref supplyMap() |
277 if ((lcap = lower[e]) != 0) { |
413 /// with such a map in which \c k is assigned to \c s, \c -k is |
278 _capacity[e] -= lcap; |
414 /// assigned to \c t and all other nodes have zero supply value. |
279 _res_cap[e] -= lcap; |
415 /// |
280 _supply[_graph.source(e)] -= lcap; |
|
281 _supply[_graph.target(e)] += lcap; |
|
282 _excess[_graph.source(e)] -= lcap; |
|
283 _excess[_graph.target(e)] += lcap; |
|
284 } |
|
285 } |
|
286 } |
|
287 /* |
|
288 /// \brief General constructor (without lower bounds). |
|
289 /// |
|
290 /// General constructor (without lower bounds). |
|
291 /// |
|
292 /// \param digraph The digraph the algorithm runs on. |
|
293 /// \param capacity The capacities (upper bounds) of the arcs. |
|
294 /// \param cost The cost (length) values of the arcs. |
|
295 /// \param supply The supply values of the nodes (signed). |
|
296 CapacityScaling( const Digraph &digraph, |
|
297 const CapacityMap &capacity, |
|
298 const CostMap &cost, |
|
299 const SupplyMap &supply ) : |
|
300 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), |
|
301 _supply(supply), _flow(NULL), _local_flow(false), |
|
302 _potential(NULL), _local_potential(false), |
|
303 _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL) |
|
304 { |
|
305 // Check the sum of supply values |
|
306 Supply sum = 0; |
|
307 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; |
|
308 _valid_supply = sum == 0; |
|
309 } |
|
310 |
|
311 /// \brief Simple constructor (with lower bounds). |
|
312 /// |
|
313 /// Simple constructor (with lower bounds). |
|
314 /// |
|
315 /// \param digraph The digraph the algorithm runs on. |
|
316 /// \param lower The lower bounds of the arcs. |
|
317 /// \param capacity The capacities (upper bounds) of the arcs. |
|
318 /// \param cost The cost (length) values of the arcs. |
|
319 /// \param s The source node. |
416 /// \param s The source node. |
320 /// \param t The target node. |
417 /// \param t The target node. |
321 /// \param flow_value The required amount of flow from node \c s |
418 /// \param k The required amount of flow from node \c s to node \c t |
322 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
419 /// (i.e. the supply of \c s and the demand of \c t). |
323 CapacityScaling( const Digraph &digraph, |
420 /// |
324 const LowerMap &lower, |
421 /// \return <tt>(*this)</tt> |
325 const CapacityMap &capacity, |
422 CapacityScaling& stSupply(const Node& s, const Node& t, Value k) { |
326 const CostMap &cost, |
423 for (int i = 0; i != _node_num; ++i) { |
327 Node s, Node t, |
424 _supply[i] = 0; |
328 Supply flow_value ) : |
425 } |
329 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost), |
426 _supply[_node_id[s]] = k; |
330 _supply(digraph, 0), _flow(NULL), _local_flow(false), |
427 _supply[_node_id[t]] = -k; |
331 _potential(NULL), _local_potential(false), |
|
332 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL) |
|
333 { |
|
334 // Remove non-zero lower bounds |
|
335 _supply[s] = _excess[s] = flow_value; |
|
336 _supply[t] = _excess[t] = -flow_value; |
|
337 typename LowerMap::Value lcap; |
|
338 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
339 if ((lcap = lower[e]) != 0) { |
|
340 _capacity[e] -= lcap; |
|
341 _res_cap[e] -= lcap; |
|
342 _supply[_graph.source(e)] -= lcap; |
|
343 _supply[_graph.target(e)] += lcap; |
|
344 _excess[_graph.source(e)] -= lcap; |
|
345 _excess[_graph.target(e)] += lcap; |
|
346 } |
|
347 } |
|
348 _valid_supply = true; |
|
349 } |
|
350 |
|
351 /// \brief Simple constructor (without lower bounds). |
|
352 /// |
|
353 /// Simple constructor (without lower bounds). |
|
354 /// |
|
355 /// \param digraph The digraph the algorithm runs on. |
|
356 /// \param capacity The capacities (upper bounds) of the arcs. |
|
357 /// \param cost The cost (length) values of the arcs. |
|
358 /// \param s The source node. |
|
359 /// \param t The target node. |
|
360 /// \param flow_value The required amount of flow from node \c s |
|
361 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
|
362 CapacityScaling( const Digraph &digraph, |
|
363 const CapacityMap &capacity, |
|
364 const CostMap &cost, |
|
365 Node s, Node t, |
|
366 Supply flow_value ) : |
|
367 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), |
|
368 _supply(digraph, 0), _flow(NULL), _local_flow(false), |
|
369 _potential(NULL), _local_potential(false), |
|
370 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL) |
|
371 { |
|
372 _supply[s] = _excess[s] = flow_value; |
|
373 _supply[t] = _excess[t] = -flow_value; |
|
374 _valid_supply = true; |
|
375 } |
|
376 */ |
|
377 /// Destructor. |
|
378 ~CapacityScaling() { |
|
379 if (_local_flow) delete _flow; |
|
380 if (_local_potential) delete _potential; |
|
381 delete _dijkstra; |
|
382 } |
|
383 |
|
384 /// \brief Set the flow map. |
|
385 /// |
|
386 /// Set the flow map. |
|
387 /// |
|
388 /// \return \c (*this) |
|
389 CapacityScaling& flowMap(FlowMap &map) { |
|
390 if (_local_flow) { |
|
391 delete _flow; |
|
392 _local_flow = false; |
|
393 } |
|
394 _flow = ↦ |
|
395 return *this; |
428 return *this; |
396 } |
429 } |
397 |
430 |
398 /// \brief Set the potential map. |
431 /// @} |
399 /// |
|
400 /// Set the potential map. |
|
401 /// |
|
402 /// \return \c (*this) |
|
403 CapacityScaling& potentialMap(PotentialMap &map) { |
|
404 if (_local_potential) { |
|
405 delete _potential; |
|
406 _local_potential = false; |
|
407 } |
|
408 _potential = ↦ |
|
409 return *this; |
|
410 } |
|
411 |
432 |
412 /// \name Execution control |
433 /// \name Execution control |
413 |
434 |
414 /// @{ |
435 /// @{ |
415 |
436 |
416 /// \brief Run the algorithm. |
437 /// \brief Run the algorithm. |
417 /// |
438 /// |
418 /// This function runs the algorithm. |
439 /// This function runs the algorithm. |
|
440 /// The paramters can be specified using functions \ref lowerMap(), |
|
441 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
|
442 /// For example, |
|
443 /// \code |
|
444 /// CapacityScaling<ListDigraph> cs(graph); |
|
445 /// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
|
446 /// .supplyMap(sup).run(); |
|
447 /// \endcode |
|
448 /// |
|
449 /// This function can be called more than once. All the parameters |
|
450 /// that have been given are kept for the next call, unless |
|
451 /// \ref reset() is called, thus only the modified parameters |
|
452 /// have to be set again. See \ref reset() for examples. |
|
453 /// However the underlying digraph must not be modified after this |
|
454 /// class have been constructed, since it copies the digraph. |
419 /// |
455 /// |
420 /// \param scaling Enable or disable capacity scaling. |
456 /// \param scaling Enable or disable capacity scaling. |
421 /// If the maximum arc capacity and/or the amount of total supply |
457 /// If the maximum upper bound and/or the amount of total supply |
422 /// is rather small, the algorithm could be slightly faster without |
458 /// is rather small, the algorithm could be slightly faster without |
423 /// scaling. |
459 /// scaling. |
424 /// |
460 /// |
425 /// \return \c true if a feasible flow can be found. |
461 /// \return \c INFEASIBLE if no feasible flow exists, |
426 bool run(bool scaling = true) { |
462 /// \n \c OPTIMAL if the problem has optimal solution |
427 return init(scaling) && start(); |
463 /// (i.e. it is feasible and bounded), and the algorithm has found |
|
464 /// optimal flow and node potentials (primal and dual solutions), |
|
465 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
|
466 /// and infinite upper bound. It means that the objective function |
|
467 /// is unbounded on that arc, however note that it could actually be |
|
468 /// bounded over the feasible flows, but this algroithm cannot handle |
|
469 /// these cases. |
|
470 /// |
|
471 /// \see ProblemType |
|
472 ProblemType run(bool scaling = true) { |
|
473 ProblemType pt = init(scaling); |
|
474 if (pt != OPTIMAL) return pt; |
|
475 return start(); |
|
476 } |
|
477 |
|
478 /// \brief Reset all the parameters that have been given before. |
|
479 /// |
|
480 /// This function resets all the paramaters that have been given |
|
481 /// before using functions \ref lowerMap(), \ref upperMap(), |
|
482 /// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
|
483 /// |
|
484 /// It is useful for multiple run() calls. If this function is not |
|
485 /// used, all the parameters given before are kept for the next |
|
486 /// \ref run() call. |
|
487 /// However the underlying digraph must not be modified after this |
|
488 /// class have been constructed, since it copies and extends the graph. |
|
489 /// |
|
490 /// For example, |
|
491 /// \code |
|
492 /// CapacityScaling<ListDigraph> cs(graph); |
|
493 /// |
|
494 /// // First run |
|
495 /// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
|
496 /// .supplyMap(sup).run(); |
|
497 /// |
|
498 /// // Run again with modified cost map (reset() is not called, |
|
499 /// // so only the cost map have to be set again) |
|
500 /// cost[e] += 100; |
|
501 /// cs.costMap(cost).run(); |
|
502 /// |
|
503 /// // Run again from scratch using reset() |
|
504 /// // (the lower bounds will be set to zero on all arcs) |
|
505 /// cs.reset(); |
|
506 /// cs.upperMap(capacity).costMap(cost) |
|
507 /// .supplyMap(sup).run(); |
|
508 /// \endcode |
|
509 /// |
|
510 /// \return <tt>(*this)</tt> |
|
511 CapacityScaling& reset() { |
|
512 for (int i = 0; i != _node_num; ++i) { |
|
513 _supply[i] = 0; |
|
514 } |
|
515 for (int j = 0; j != _res_arc_num; ++j) { |
|
516 _lower[j] = 0; |
|
517 _upper[j] = INF; |
|
518 _cost[j] = _forward[j] ? 1 : -1; |
|
519 } |
|
520 _have_lower = false; |
|
521 return *this; |
428 } |
522 } |
429 |
523 |
430 /// @} |
524 /// @} |
431 |
525 |
432 /// \name Query Functions |
526 /// \name Query Functions |
433 /// The results of the algorithm can be obtained using these |
527 /// The results of the algorithm can be obtained using these |
434 /// functions.\n |
528 /// functions.\n |
435 /// \ref lemon::CapacityScaling::run() "run()" must be called before |
529 /// The \ref run() function must be called before using them. |
436 /// using them. |
|
437 |
530 |
438 /// @{ |
531 /// @{ |
439 |
532 |
440 /// \brief Return a const reference to the arc map storing the |
533 /// \brief Return the total cost of the found flow. |
441 /// found flow. |
534 /// |
442 /// |
535 /// This function returns the total cost of the found flow. |
443 /// Return a const reference to the arc map storing the found flow. |
536 /// Its complexity is O(e). |
|
537 /// |
|
538 /// \note The return type of the function can be specified as a |
|
539 /// template parameter. For example, |
|
540 /// \code |
|
541 /// cs.totalCost<double>(); |
|
542 /// \endcode |
|
543 /// It is useful if the total cost cannot be stored in the \c Cost |
|
544 /// type of the algorithm, which is the default return type of the |
|
545 /// function. |
444 /// |
546 /// |
445 /// \pre \ref run() must be called before using this function. |
547 /// \pre \ref run() must be called before using this function. |
446 const FlowMap& flowMap() const { |
548 template <typename Number> |
447 return *_flow; |
549 Number totalCost() const { |
448 } |
550 Number c = 0; |
449 |
551 for (ArcIt a(_graph); a != INVALID; ++a) { |
450 /// \brief Return a const reference to the node map storing the |
552 int i = _arc_idb[a]; |
451 /// found potentials (the dual solution). |
553 c += static_cast<Number>(_res_cap[i]) * |
452 /// |
554 (-static_cast<Number>(_cost[i])); |
453 /// Return a const reference to the node map storing the found |
555 } |
454 /// potentials (the dual solution). |
556 return c; |
|
557 } |
|
558 |
|
559 #ifndef DOXYGEN |
|
560 Cost totalCost() const { |
|
561 return totalCost<Cost>(); |
|
562 } |
|
563 #endif |
|
564 |
|
565 /// \brief Return the flow on the given arc. |
|
566 /// |
|
567 /// This function returns the flow on the given arc. |
455 /// |
568 /// |
456 /// \pre \ref run() must be called before using this function. |
569 /// \pre \ref run() must be called before using this function. |
457 const PotentialMap& potentialMap() const { |
570 Value flow(const Arc& a) const { |
458 return *_potential; |
571 return _res_cap[_arc_idb[a]]; |
459 } |
572 } |
460 |
573 |
461 /// \brief Return the flow on the given arc. |
574 /// \brief Return the flow map (the primal solution). |
462 /// |
575 /// |
463 /// Return the flow on the given arc. |
576 /// This function copies the flow value on each arc into the given |
|
577 /// map. The \c Value type of the algorithm must be convertible to |
|
578 /// the \c Value type of the map. |
464 /// |
579 /// |
465 /// \pre \ref run() must be called before using this function. |
580 /// \pre \ref run() must be called before using this function. |
466 Capacity flow(const Arc& arc) const { |
581 template <typename FlowMap> |
467 return (*_flow)[arc]; |
582 void flowMap(FlowMap &map) const { |
468 } |
583 for (ArcIt a(_graph); a != INVALID; ++a) { |
469 |
584 map.set(a, _res_cap[_arc_idb[a]]); |
470 /// \brief Return the potential of the given node. |
585 } |
471 /// |
586 } |
472 /// Return the potential of the given node. |
587 |
|
588 /// \brief Return the potential (dual value) of the given node. |
|
589 /// |
|
590 /// This function returns the potential (dual value) of the |
|
591 /// given node. |
473 /// |
592 /// |
474 /// \pre \ref run() must be called before using this function. |
593 /// \pre \ref run() must be called before using this function. |
475 Cost potential(const Node& node) const { |
594 Cost potential(const Node& n) const { |
476 return (*_potential)[node]; |
595 return _pi[_node_id[n]]; |
477 } |
596 } |
478 |
597 |
479 /// \brief Return the total cost of the found flow. |
598 /// \brief Return the potential map (the dual solution). |
480 /// |
599 /// |
481 /// Return the total cost of the found flow. The complexity of the |
600 /// This function copies the potential (dual value) of each node |
482 /// function is \f$ O(e) \f$. |
601 /// into the given map. |
|
602 /// The \c Cost type of the algorithm must be convertible to the |
|
603 /// \c Value type of the map. |
483 /// |
604 /// |
484 /// \pre \ref run() must be called before using this function. |
605 /// \pre \ref run() must be called before using this function. |
485 Cost totalCost() const { |
606 template <typename PotentialMap> |
486 Cost c = 0; |
607 void potentialMap(PotentialMap &map) const { |
487 for (ArcIt e(_graph); e != INVALID; ++e) |
608 for (NodeIt n(_graph); n != INVALID; ++n) { |
488 c += (*_flow)[e] * _cost[e]; |
609 map.set(n, _pi[_node_id[n]]); |
489 return c; |
610 } |
490 } |
611 } |
491 |
612 |
492 /// @} |
613 /// @} |
493 |
614 |
494 private: |
615 private: |
495 |
616 |
496 /// Initialize the algorithm. |
617 // Initialize the algorithm |
497 bool init(bool scaling) { |
618 ProblemType init(bool scaling) { |
498 if (!_valid_supply) return false; |
619 if (_node_num == 0) return INFEASIBLE; |
499 |
620 |
500 // Initializing maps |
621 // Check the sum of supply values |
501 if (!_flow) { |
622 _sum_supply = 0; |
502 _flow = new FlowMap(_graph); |
623 for (int i = 0; i != _root; ++i) { |
503 _local_flow = true; |
624 _sum_supply += _supply[i]; |
504 } |
625 } |
505 if (!_potential) { |
626 if (_sum_supply > 0) return INFEASIBLE; |
506 _potential = new PotentialMap(_graph); |
627 |
507 _local_potential = true; |
628 // Initialize maps |
508 } |
629 for (int i = 0; i != _root; ++i) { |
509 for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0; |
630 _pi[i] = 0; |
510 for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0; |
631 _excess[i] = _supply[i]; |
511 |
632 } |
512 _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost, |
633 |
513 _excess, *_potential, _pred ); |
634 // Remove non-zero lower bounds |
514 |
635 if (_have_lower) { |
515 // Initializing delta value |
636 for (int i = 0; i != _root; ++i) { |
|
637 for (int j = _first_out[i]; j != _first_out[i+1]; ++j) { |
|
638 if (_forward[j]) { |
|
639 Value c = _lower[j]; |
|
640 if (c >= 0) { |
|
641 _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF; |
|
642 } else { |
|
643 _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF; |
|
644 } |
|
645 _excess[i] -= c; |
|
646 _excess[_target[j]] += c; |
|
647 } else { |
|
648 _res_cap[j] = 0; |
|
649 } |
|
650 } |
|
651 } |
|
652 } else { |
|
653 for (int j = 0; j != _res_arc_num; ++j) { |
|
654 _res_cap[j] = _forward[j] ? _upper[j] : 0; |
|
655 } |
|
656 } |
|
657 |
|
658 // Handle negative costs |
|
659 for (int u = 0; u != _root; ++u) { |
|
660 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { |
|
661 Value rc = _res_cap[a]; |
|
662 if (_cost[a] < 0 && rc > 0) { |
|
663 if (rc == INF) return UNBOUNDED; |
|
664 _excess[u] -= rc; |
|
665 _excess[_target[a]] += rc; |
|
666 _res_cap[a] = 0; |
|
667 _res_cap[_reverse[a]] += rc; |
|
668 } |
|
669 } |
|
670 } |
|
671 |
|
672 // Handle GEQ supply type |
|
673 if (_sum_supply < 0) { |
|
674 _pi[_root] = 0; |
|
675 _excess[_root] = -_sum_supply; |
|
676 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
|
677 int u = _target[a]; |
|
678 if (_excess[u] < 0) { |
|
679 _res_cap[a] = -_excess[u] + 1; |
|
680 } else { |
|
681 _res_cap[a] = 1; |
|
682 } |
|
683 _res_cap[_reverse[a]] = 0; |
|
684 _cost[a] = 0; |
|
685 _cost[_reverse[a]] = 0; |
|
686 } |
|
687 } else { |
|
688 _pi[_root] = 0; |
|
689 _excess[_root] = 0; |
|
690 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
|
691 _res_cap[a] = 1; |
|
692 _res_cap[_reverse[a]] = 0; |
|
693 _cost[a] = 0; |
|
694 _cost[_reverse[a]] = 0; |
|
695 } |
|
696 } |
|
697 |
|
698 // Initialize delta value |
516 if (scaling) { |
699 if (scaling) { |
517 // With scaling |
700 // With scaling |
518 Supply max_sup = 0, max_dem = 0; |
701 Value max_sup = 0, max_dem = 0; |
519 for (NodeIt n(_graph); n != INVALID; ++n) { |
702 for (int i = 0; i != _node_num; ++i) { |
520 if ( _supply[n] > max_sup) max_sup = _supply[n]; |
703 if ( _excess[i] > max_sup) max_sup = _excess[i]; |
521 if (-_supply[n] > max_dem) max_dem = -_supply[n]; |
704 if (-_excess[i] > max_dem) max_dem = -_excess[i]; |
522 } |
705 } |
523 Capacity max_cap = 0; |
706 Value max_cap = 0; |
524 for (ArcIt e(_graph); e != INVALID; ++e) { |
707 for (int j = 0; j != _res_arc_num; ++j) { |
525 if (_capacity[e] > max_cap) max_cap = _capacity[e]; |
708 if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; |
526 } |
709 } |
527 max_sup = std::min(std::min(max_sup, max_dem), max_cap); |
710 max_sup = std::min(std::min(max_sup, max_dem), max_cap); |
528 _phase_num = 0; |
711 _phase_num = 0; |
529 for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) |
712 for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) |
530 ++_phase_num; |
713 ++_phase_num; |
531 } else { |
714 } else { |
532 // Without scaling |
715 // Without scaling |
533 _delta = 1; |
716 _delta = 1; |
534 } |
717 } |
535 |
718 |
536 return true; |
719 return OPTIMAL; |
537 } |
720 } |
538 |
721 |
539 bool start() { |
722 ProblemType start() { |
|
723 // Execute the algorithm |
|
724 ProblemType pt; |
540 if (_delta > 1) |
725 if (_delta > 1) |
541 return startWithScaling(); |
726 pt = startWithScaling(); |
542 else |
727 else |
543 return startWithoutScaling(); |
728 pt = startWithoutScaling(); |
544 } |
729 |
545 |
730 // Handle non-zero lower bounds |
546 /// Execute the capacity scaling algorithm. |
731 if (_have_lower) { |
547 bool startWithScaling() { |
732 for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) { |
548 // Processing capacity scaling phases |
733 if (!_forward[j]) _res_cap[j] += _lower[j]; |
549 Node s, t; |
734 } |
|
735 } |
|
736 |
|
737 // Shift potentials if necessary |
|
738 Cost pr = _pi[_root]; |
|
739 if (_sum_supply < 0 || pr > 0) { |
|
740 for (int i = 0; i != _node_num; ++i) { |
|
741 _pi[i] -= pr; |
|
742 } |
|
743 } |
|
744 |
|
745 return pt; |
|
746 } |
|
747 |
|
748 // Execute the capacity scaling algorithm |
|
749 ProblemType startWithScaling() { |
|
750 // Process capacity scaling phases |
|
751 int s, t; |
550 int phase_cnt = 0; |
752 int phase_cnt = 0; |
551 int factor = 4; |
753 int factor = 4; |
|
754 ResidualDijkstra _dijkstra(*this); |
552 while (true) { |
755 while (true) { |
553 // Saturating all arcs not satisfying the optimality condition |
756 // Saturate all arcs not satisfying the optimality condition |
554 for (ArcIt e(_graph); e != INVALID; ++e) { |
757 for (int u = 0; u != _node_num; ++u) { |
555 Node u = _graph.source(e), v = _graph.target(e); |
758 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { |
556 Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v]; |
759 int v = _target[a]; |
557 if (c < 0 && _res_cap[e] >= _delta) { |
760 Cost c = _cost[a] + _pi[u] - _pi[v]; |
558 _excess[u] -= _res_cap[e]; |
761 Value rc = _res_cap[a]; |
559 _excess[v] += _res_cap[e]; |
762 if (c < 0 && rc >= _delta) { |
560 (*_flow)[e] = _capacity[e]; |
763 _excess[u] -= rc; |
561 _res_cap[e] = 0; |
764 _excess[v] += rc; |
562 } |
765 _res_cap[a] = 0; |
563 else if (c > 0 && (*_flow)[e] >= _delta) { |
766 _res_cap[_reverse[a]] += rc; |
564 _excess[u] += (*_flow)[e]; |
767 } |
565 _excess[v] -= (*_flow)[e]; |
768 } |
566 (*_flow)[e] = 0; |
769 } |
567 _res_cap[e] = _capacity[e]; |
770 |
568 } |
771 // Find excess nodes and deficit nodes |
569 } |
|
570 |
|
571 // Finding excess nodes and deficit nodes |
|
572 _excess_nodes.clear(); |
772 _excess_nodes.clear(); |
573 _deficit_nodes.clear(); |
773 _deficit_nodes.clear(); |
574 for (NodeIt n(_graph); n != INVALID; ++n) { |
774 for (int u = 0; u != _node_num; ++u) { |
575 if (_excess[n] >= _delta) _excess_nodes.push_back(n); |
775 if (_excess[u] >= _delta) _excess_nodes.push_back(u); |
576 if (_excess[n] <= -_delta) _deficit_nodes.push_back(n); |
776 if (_excess[u] <= -_delta) _deficit_nodes.push_back(u); |
577 } |
777 } |
578 int next_node = 0, next_def_node = 0; |
778 int next_node = 0, next_def_node = 0; |
579 |
779 |
580 // Finding augmenting shortest paths |
780 // Find augmenting shortest paths |
581 while (next_node < int(_excess_nodes.size())) { |
781 while (next_node < int(_excess_nodes.size())) { |
582 // Checking deficit nodes |
782 // Check deficit nodes |
583 if (_delta > 1) { |
783 if (_delta > 1) { |
584 bool delta_deficit = false; |
784 bool delta_deficit = false; |
585 for ( ; next_def_node < int(_deficit_nodes.size()); |
785 for ( ; next_def_node < int(_deficit_nodes.size()); |
586 ++next_def_node ) { |
786 ++next_def_node ) { |
587 if (_excess[_deficit_nodes[next_def_node]] <= -_delta) { |
787 if (_excess[_deficit_nodes[next_def_node]] <= -_delta) { |