23 /// |
23 /// |
24 /// \file |
24 /// \file |
25 /// \brief The capacity scaling algorithm for finding a minimum cost |
25 /// \brief The capacity scaling algorithm for finding a minimum cost |
26 /// flow. |
26 /// flow. |
27 |
27 |
|
28 #include <lemon/graph_adaptor.h> |
|
29 #include <lemon/bin_heap.h> |
28 #include <vector> |
30 #include <vector> |
29 #include <lemon/graph_adaptor.h> |
|
30 #include <lemon/dijkstra.h> |
|
31 |
|
32 #define WITH_SCALING |
|
33 |
|
34 #ifdef WITH_SCALING |
|
35 #define SCALING_FACTOR 2 |
|
36 #endif |
|
37 |
|
38 //#define _DEBUG_ITER_ |
|
39 |
31 |
40 namespace lemon { |
32 namespace lemon { |
41 |
33 |
42 /// \addtogroup min_cost_flow |
34 /// \addtogroup min_cost_flow |
43 /// @{ |
35 /// @{ |
44 |
36 |
45 /// \brief Implementation of the capacity scaling version of the |
37 /// \brief Implementation of the capacity scaling version of the |
46 /// successive shortest path algorithm for finding a minimum cost |
38 /// successive shortest path algorithm for finding a minimum cost |
47 /// flow. |
39 /// flow. |
48 /// |
40 /// |
49 /// \ref lemon::CapacityScaling "CapacityScaling" implements the |
41 /// \ref CapacityScaling implements the capacity scaling version |
50 /// capacity scaling version of the successive shortest path |
42 /// of the successive shortest path algorithm for finding a minimum |
51 /// algorithm for finding a minimum cost flow. |
43 /// cost flow. |
52 /// |
44 /// |
53 /// \param Graph The directed graph type the algorithm runs on. |
45 /// \param Graph The directed graph type the algorithm runs on. |
54 /// \param LowerMap The type of the lower bound map. |
46 /// \param LowerMap The type of the lower bound map. |
55 /// \param CapacityMap The type of the capacity (upper bound) map. |
47 /// \param CapacityMap The type of the capacity (upper bound) map. |
56 /// \param CostMap The type of the cost (length) map. |
48 /// \param CostMap The type of the cost (length) map. |
57 /// \param SupplyMap The type of the supply map. |
49 /// \param SupplyMap The type of the supply map. |
58 /// |
50 /// |
59 /// \warning |
51 /// \warning |
60 /// - Edge capacities and costs should be nonnegative integers. |
52 /// - Edge capacities and costs should be nonnegative integers. |
61 /// However \c CostMap::Value should be signed type. |
53 /// However \c CostMap::Value should be signed type. |
62 /// - Supply values should be signed integers. |
54 /// - Supply values should be signed integers. |
63 /// - \c LowerMap::Value must be convertible to |
55 /// - \c LowerMap::Value must be convertible to |
64 /// \c CapacityMap::Value and \c CapacityMap::Value must be |
56 /// \c CapacityMap::Value and \c CapacityMap::Value must be |
65 /// convertible to \c SupplyMap::Value. |
57 /// convertible to \c SupplyMap::Value. |
66 /// |
58 /// |
67 /// \author Peter Kovacs |
59 /// \author Peter Kovacs |
68 |
60 |
69 template < typename Graph, |
61 template < typename Graph, |
70 typename LowerMap = typename Graph::template EdgeMap<int>, |
62 typename LowerMap = typename Graph::template EdgeMap<int>, |
71 typename CapacityMap = LowerMap, |
63 typename CapacityMap = LowerMap, |
72 typename CostMap = typename Graph::template EdgeMap<int>, |
64 typename CostMap = typename Graph::template EdgeMap<int>, |
73 typename SupplyMap = typename Graph::template NodeMap |
65 typename SupplyMap = typename Graph::template NodeMap |
74 <typename CapacityMap::Value> > |
66 <typename CapacityMap::Value> > |
75 class CapacityScaling |
67 class CapacityScaling |
76 { |
68 { |
77 typedef typename Graph::Node Node; |
69 typedef typename Graph::Node Node; |
78 typedef typename Graph::NodeIt NodeIt; |
70 typedef typename Graph::NodeIt NodeIt; |
79 typedef typename Graph::Edge Edge; |
71 typedef typename Graph::Edge Edge; |
85 typedef typename CapacityMap::Value Capacity; |
77 typedef typename CapacityMap::Value Capacity; |
86 typedef typename CostMap::Value Cost; |
78 typedef typename CostMap::Value Cost; |
87 typedef typename SupplyMap::Value Supply; |
79 typedef typename SupplyMap::Value Supply; |
88 typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap; |
80 typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap; |
89 typedef typename Graph::template NodeMap<Supply> SupplyRefMap; |
81 typedef typename Graph::template NodeMap<Supply> SupplyRefMap; |
90 |
82 typedef typename Graph::template NodeMap<Edge> PredMap; |
91 typedef ResGraphAdaptor< const Graph, Capacity, |
|
92 CapacityRefMap, CapacityRefMap > ResGraph; |
|
93 typedef typename ResGraph::Node ResNode; |
|
94 typedef typename ResGraph::NodeIt ResNodeIt; |
|
95 typedef typename ResGraph::Edge ResEdge; |
|
96 typedef typename ResGraph::EdgeIt ResEdgeIt; |
|
97 |
83 |
98 public: |
84 public: |
|
85 |
|
86 /// \brief Type to enable or disable capacity scaling. |
|
87 enum ScalingEnum { |
|
88 WITH_SCALING = 0, |
|
89 WITHOUT_SCALING = -1 |
|
90 }; |
99 |
91 |
100 /// \brief The type of the flow map. |
92 /// \brief The type of the flow map. |
101 typedef CapacityRefMap FlowMap; |
93 typedef CapacityRefMap FlowMap; |
102 /// \brief The type of the potential map. |
94 /// \brief The type of the potential map. |
103 typedef typename Graph::template NodeMap<Cost> PotentialMap; |
95 typedef typename Graph::template NodeMap<Cost> PotentialMap; |
104 |
96 |
105 protected: |
97 protected: |
106 |
98 |
107 /// \brief Map adaptor class for handling reduced edge costs. |
99 /// \brief Special implementation of the \ref Dijkstra algorithm |
108 class ReducedCostMap : public MapBase<ResEdge, Cost> |
100 /// for finding shortest paths in the residual network of the graph |
|
101 /// with respect to the reduced edge costs and modifying the |
|
102 /// node potentials according to the distance of the nodes. |
|
103 class ResidualDijkstra |
109 { |
104 { |
110 private: |
105 typedef typename Graph::template NodeMap<Cost> CostNodeMap; |
111 |
106 typedef typename Graph::template NodeMap<Edge> PredMap; |
112 const ResGraph &gr; |
107 |
113 const CostMap &cost_map; |
108 typedef typename Graph::template NodeMap<int> HeapCrossRef; |
114 const PotentialMap &pot_map; |
109 typedef BinHeap<Cost, HeapCrossRef> Heap; |
|
110 |
|
111 protected: |
|
112 |
|
113 /// \brief The directed graph the algorithm runs on. |
|
114 const Graph &graph; |
|
115 |
|
116 /// \brief The flow map. |
|
117 const FlowMap &flow; |
|
118 /// \brief The residual capacity map. |
|
119 const CapacityRefMap &res_cap; |
|
120 /// \brief The cost map. |
|
121 const CostMap &cost; |
|
122 /// \brief The excess map. |
|
123 const SupplyRefMap &excess; |
|
124 |
|
125 /// \brief The potential map. |
|
126 PotentialMap &potential; |
|
127 |
|
128 /// \brief The distance map. |
|
129 CostNodeMap dist; |
|
130 /// \brief The map of predecessors edges. |
|
131 PredMap &pred; |
|
132 /// \brief The processed (i.e. permanently labeled) nodes. |
|
133 std::vector<Node> proc_nodes; |
115 |
134 |
116 public: |
135 public: |
117 |
136 |
118 ReducedCostMap( const ResGraph &_gr, |
137 /// \brief The constructor of the class. |
119 const CostMap &_cost, |
138 ResidualDijkstra( const Graph &_graph, |
120 const PotentialMap &_pot ) : |
139 const FlowMap &_flow, |
121 gr(_gr), cost_map(_cost), pot_map(_pot) {} |
140 const CapacityRefMap &_res_cap, |
122 |
141 const CostMap &_cost, |
123 Cost operator[](const ResEdge &e) const { |
142 const SupplyMap &_excess, |
124 return ResGraph::forward(e) ? |
143 PotentialMap &_potential, |
125 cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] : |
144 PredMap &_pred ) : |
126 -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)]; |
145 graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost), |
127 } |
146 excess(_excess), potential(_potential), dist(_graph), |
128 |
147 pred(_pred) |
129 }; //class ReducedCostMap |
148 {} |
130 |
149 |
131 /// \brief Map class for the \ref lemon::Dijkstra "Dijkstra" |
150 /// \brief Runs the algorithm from the given source node. |
132 /// algorithm to update node potentials. |
151 Node run(Node s, Capacity delta) { |
133 class PotentialUpdateMap : public MapBase<ResNode, Cost> |
152 HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP); |
134 { |
153 Heap heap(heap_cross_ref); |
135 private: |
154 heap.push(s, 0); |
136 |
155 pred[s] = INVALID; |
137 PotentialMap *pot; |
156 proc_nodes.clear(); |
138 typedef std::pair<ResNode, Cost> Pair; |
157 |
139 std::vector<Pair> data; |
158 // Processing nodes |
140 |
159 while (!heap.empty() && excess[heap.top()] > -delta) { |
141 public: |
160 Node u = heap.top(), v; |
142 |
161 Cost d = heap.prio() - potential[u], nd; |
143 void potentialMap(PotentialMap &_pot) { |
162 dist[u] = heap.prio(); |
144 pot = &_pot; |
163 heap.pop(); |
145 } |
164 proc_nodes.push_back(u); |
146 |
165 |
147 void init() { |
166 // Traversing outgoing edges |
148 data.clear(); |
167 for (OutEdgeIt e(graph, u); e != INVALID; ++e) { |
149 } |
168 if (res_cap[e] >= delta) { |
150 |
169 v = graph.target(e); |
151 void set(const ResNode &n, const Cost &v) { |
170 switch(heap.state(v)) { |
152 data.push_back(Pair(n, v)); |
171 case Heap::PRE_HEAP: |
153 } |
172 heap.push(v, d + cost[e] + potential[v]); |
154 |
173 pred[v] = e; |
155 void update() { |
174 break; |
156 Cost val = data[data.size()-1].second; |
175 case Heap::IN_HEAP: |
157 for (int i = 0; i < data.size()-1; ++i) |
176 nd = d + cost[e] + potential[v]; |
158 (*pot)[data[i].first] += val - data[i].second; |
177 if (nd < heap[v]) { |
159 } |
178 heap.decrease(v, nd); |
160 |
179 pred[v] = e; |
161 }; //class PotentialUpdateMap |
180 } |
162 |
181 break; |
163 #ifdef WITH_SCALING |
182 case Heap::POST_HEAP: |
164 /// \brief Map adaptor class for identifing deficit nodes. |
183 break; |
165 class DeficitBoolMap : public MapBase<ResNode, bool> |
184 } |
166 { |
185 } |
167 private: |
186 } |
168 |
187 |
169 const SupplyRefMap &imb; |
188 // Traversing incoming edges |
170 const Capacity δ |
189 for (InEdgeIt e(graph, u); e != INVALID; ++e) { |
171 |
190 if (flow[e] >= delta) { |
172 public: |
191 v = graph.source(e); |
173 |
192 switch(heap.state(v)) { |
174 DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) : |
193 case Heap::PRE_HEAP: |
175 imb(_imb), delta(_delta) {} |
194 heap.push(v, d - cost[e] + potential[v]); |
176 |
195 pred[v] = e; |
177 bool operator[](const ResNode &n) const { |
196 break; |
178 return imb[n] <= -delta; |
197 case Heap::IN_HEAP: |
179 } |
198 nd = d - cost[e] + potential[v]; |
180 |
199 if (nd < heap[v]) { |
181 }; //class DeficitBoolMap |
200 heap.decrease(v, nd); |
182 |
201 pred[v] = e; |
183 /// \brief Map adaptor class for filtering edges with at least |
202 } |
184 /// \c delta residual capacity. |
203 break; |
185 class DeltaFilterMap : public MapBase<ResEdge, bool> |
204 case Heap::POST_HEAP: |
186 { |
205 break; |
187 private: |
206 } |
188 |
207 } |
189 const ResGraph &gr; |
208 } |
190 const Capacity δ |
209 } |
191 |
210 if (heap.empty()) return INVALID; |
192 public: |
211 |
193 |
212 // Updating potentials of processed nodes |
194 DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) : |
213 Node t = heap.top(); |
195 gr(_gr), delta(_delta) {} |
214 Cost dt = heap.prio(); |
196 |
215 for (int i = 0; i < proc_nodes.size(); ++i) |
197 bool operator[](const ResEdge &e) const { |
216 potential[proc_nodes[i]] -= dist[proc_nodes[i]] - dt; |
198 return gr.rescap(e) >= delta; |
217 |
199 } |
218 return t; |
200 |
219 } |
201 }; //class DeltaFilterMap |
220 |
202 |
221 }; //class ResidualDijkstra |
203 typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap> |
|
204 DeltaResGraph; |
|
205 |
|
206 /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class. |
|
207 class ResDijkstraTraits : |
|
208 public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap> |
|
209 { |
|
210 public: |
|
211 |
|
212 typedef PotentialUpdateMap DistMap; |
|
213 |
|
214 static DistMap *createDistMap(const DeltaResGraph&) { |
|
215 return new DistMap(); |
|
216 } |
|
217 |
|
218 }; //class ResDijkstraTraits |
|
219 |
|
220 #else //WITHOUT_CAPACITY_SCALING |
|
221 /// \brief Map adaptor class for identifing deficit nodes. |
|
222 class DeficitBoolMap : public MapBase<ResNode, bool> |
|
223 { |
|
224 private: |
|
225 |
|
226 const SupplyRefMap &imb; |
|
227 |
|
228 public: |
|
229 |
|
230 DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {} |
|
231 |
|
232 bool operator[](const ResNode &n) const { |
|
233 return imb[n] < 0; |
|
234 } |
|
235 |
|
236 }; //class DeficitBoolMap |
|
237 |
|
238 /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class. |
|
239 class ResDijkstraTraits : |
|
240 public DijkstraDefaultTraits<ResGraph, ReducedCostMap> |
|
241 { |
|
242 public: |
|
243 |
|
244 typedef PotentialUpdateMap DistMap; |
|
245 |
|
246 static DistMap *createDistMap(const ResGraph&) { |
|
247 return new DistMap(); |
|
248 } |
|
249 |
|
250 }; //class ResDijkstraTraits |
|
251 #endif |
|
252 |
222 |
253 protected: |
223 protected: |
254 |
224 |
255 /// \brief The directed graph the algorithm runs on. |
225 /// \brief The directed graph the algorithm runs on. |
256 const Graph &graph; |
226 const Graph &graph; |
267 |
237 |
268 /// \brief The edge map of the current flow. |
238 /// \brief The edge map of the current flow. |
269 FlowMap flow; |
239 FlowMap flow; |
270 /// \brief The potential node map. |
240 /// \brief The potential node map. |
271 PotentialMap potential; |
241 PotentialMap potential; |
272 /// \brief The residual graph. |
242 |
273 ResGraph res_graph; |
243 /// \brief The residual capacity map. |
274 /// \brief The reduced cost map. |
244 CapacityRefMap res_cap; |
275 ReducedCostMap red_cost; |
245 /// \brief The excess map. |
276 |
246 SupplyRefMap excess; |
277 /// \brief The imbalance map. |
247 /// \brief The excess nodes (i.e. nodes with positive excess). |
278 SupplyRefMap imbalance; |
248 std::vector<Node> excess_nodes; |
279 /// \brief The excess nodes. |
|
280 std::vector<ResNode> excess_nodes; |
|
281 /// \brief The index of the next excess node. |
249 /// \brief The index of the next excess node. |
282 int next_node; |
250 int next_node; |
283 |
251 |
284 #ifdef WITH_SCALING |
252 /// \brief The scaling status (enabled or disabled). |
285 typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits> |
253 ScalingEnum scaling; |
286 ResDijkstra; |
|
287 /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding |
|
288 /// shortest paths in the residual graph with respect to the |
|
289 /// reduced edge costs. |
|
290 ResDijkstra dijkstra; |
|
291 |
|
292 /// \brief The delta parameter used for capacity scaling. |
254 /// \brief The delta parameter used for capacity scaling. |
293 Capacity delta; |
255 Capacity delta; |
294 /// \brief Edge filter map for filtering edges with at least |
256 /// \brief The maximum number of phases. |
295 /// \c delta residual capacity. |
257 Capacity phase_num; |
296 DeltaFilterMap delta_filter; |
|
297 /// \brief The delta residual graph (i.e. the subgraph of the |
|
298 /// residual graph consisting of edges with at least \c delta |
|
299 /// residual capacity). |
|
300 DeltaResGraph dres_graph; |
|
301 /// \brief Map for identifing deficit nodes. |
|
302 DeficitBoolMap delta_deficit; |
|
303 |
|
304 /// \brief The deficit nodes. |
258 /// \brief The deficit nodes. |
305 std::vector<ResNode> deficit_nodes; |
259 std::vector<Node> deficit_nodes; |
306 |
260 |
307 #else //WITHOUT_CAPACITY_SCALING |
261 /// \brief Implementation of the \ref Dijkstra algorithm for |
308 typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits> |
262 /// finding augmenting shortest paths in the residual network. |
309 ResDijkstra; |
263 ResidualDijkstra dijkstra; |
310 /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding |
264 /// \brief The map of predecessors edges. |
311 /// shortest paths in the residual graph with respect to the |
265 PredMap pred; |
312 /// reduced edge costs. |
|
313 ResDijkstra dijkstra; |
|
314 /// \brief Map for identifing deficit nodes. |
|
315 DeficitBoolMap has_deficit; |
|
316 #endif |
|
317 |
|
318 /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class. |
|
319 typename ResDijkstra::PredMap pred; |
|
320 /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to |
|
321 /// update node potentials. |
|
322 PotentialUpdateMap updater; |
|
323 |
266 |
324 public : |
267 public : |
325 |
268 |
326 /// \brief General constructor of the class (with lower bounds). |
269 /// \brief General constructor of the class (with lower bounds). |
327 /// |
270 /// |
331 /// \param _lower The lower bounds of the edges. |
274 /// \param _lower The lower bounds of the edges. |
332 /// \param _capacity The capacities (upper bounds) of the edges. |
275 /// \param _capacity The capacities (upper bounds) of the edges. |
333 /// \param _cost The cost (length) values of the edges. |
276 /// \param _cost The cost (length) values of the edges. |
334 /// \param _supply The supply values of the nodes (signed). |
277 /// \param _supply The supply values of the nodes (signed). |
335 CapacityScaling( const Graph &_graph, |
278 CapacityScaling( const Graph &_graph, |
336 const LowerMap &_lower, |
279 const LowerMap &_lower, |
337 const CapacityMap &_capacity, |
280 const CapacityMap &_capacity, |
338 const CostMap &_cost, |
281 const CostMap &_cost, |
339 const SupplyMap &_supply ) : |
282 const SupplyMap &_supply ) : |
340 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost), |
283 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost), |
341 supply(_graph), flow(_graph, 0), potential(_graph, 0), |
284 supply(_graph), flow(_graph, 0), potential(_graph, 0), |
342 res_graph(_graph, capacity, flow), |
285 res_cap(_graph), excess(_graph), pred(_graph), |
343 red_cost(res_graph, cost, potential), imbalance(_graph), |
286 dijkstra(graph, flow, res_cap, cost, excess, potential, pred) |
344 #ifdef WITH_SCALING |
|
345 delta(0), delta_filter(res_graph, delta), |
|
346 dres_graph(res_graph, delta_filter), |
|
347 dijkstra(dres_graph, red_cost), pred(dres_graph), |
|
348 delta_deficit(imbalance, delta) |
|
349 #else //WITHOUT_CAPACITY_SCALING |
|
350 dijkstra(res_graph, red_cost), pred(res_graph), |
|
351 has_deficit(imbalance) |
|
352 #endif |
|
353 { |
287 { |
354 // Removing nonzero lower bounds |
288 // Removing nonzero lower bounds |
355 capacity = subMap(_capacity, _lower); |
289 capacity = subMap(_capacity, _lower); |
|
290 res_cap = capacity; |
356 Supply sum = 0; |
291 Supply sum = 0; |
357 for (NodeIt n(graph); n != INVALID; ++n) { |
292 for (NodeIt n(graph); n != INVALID; ++n) { |
358 Supply s = _supply[n]; |
293 Supply s = _supply[n]; |
359 for (InEdgeIt e(graph, n); e != INVALID; ++e) |
294 for (InEdgeIt e(graph, n); e != INVALID; ++e) |
360 s += _lower[e]; |
295 s += _lower[e]; |
361 for (OutEdgeIt e(graph, n); e != INVALID; ++e) |
296 for (OutEdgeIt e(graph, n); e != INVALID; ++e) |
362 s -= _lower[e]; |
297 s -= _lower[e]; |
363 supply[n] = s; |
298 supply[n] = s; |
364 sum += s; |
299 sum += s; |
365 } |
300 } |
366 valid_supply = sum == 0; |
301 valid_supply = sum == 0; |
367 } |
302 } |
368 |
303 |
369 /// \brief General constructor of the class (without lower bounds). |
304 /// \brief General constructor of the class (without lower bounds). |
373 /// \param _graph The directed graph the algorithm runs on. |
308 /// \param _graph The directed graph the algorithm runs on. |
374 /// \param _capacity The capacities (upper bounds) of the edges. |
309 /// \param _capacity The capacities (upper bounds) of the edges. |
375 /// \param _cost The cost (length) values of the edges. |
310 /// \param _cost The cost (length) values of the edges. |
376 /// \param _supply The supply values of the nodes (signed). |
311 /// \param _supply The supply values of the nodes (signed). |
377 CapacityScaling( const Graph &_graph, |
312 CapacityScaling( const Graph &_graph, |
378 const CapacityMap &_capacity, |
313 const CapacityMap &_capacity, |
379 const CostMap &_cost, |
314 const CostMap &_cost, |
380 const SupplyMap &_supply ) : |
315 const SupplyMap &_supply ) : |
381 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost), |
316 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost), |
382 supply(_supply), flow(_graph, 0), potential(_graph, 0), |
317 supply(_supply), flow(_graph, 0), potential(_graph, 0), |
383 res_graph(_graph, capacity, flow), |
318 res_cap(_capacity), excess(_graph), pred(_graph), |
384 red_cost(res_graph, cost, potential), imbalance(_graph), |
319 dijkstra(graph, flow, res_cap, cost, excess, potential) |
385 #ifdef WITH_SCALING |
|
386 delta(0), delta_filter(res_graph, delta), |
|
387 dres_graph(res_graph, delta_filter), |
|
388 dijkstra(dres_graph, red_cost), pred(dres_graph), |
|
389 delta_deficit(imbalance, delta) |
|
390 #else //WITHOUT_CAPACITY_SCALING |
|
391 dijkstra(res_graph, red_cost), pred(res_graph), |
|
392 has_deficit(imbalance) |
|
393 #endif |
|
394 { |
320 { |
395 // Checking the sum of supply values |
321 // Checking the sum of supply values |
396 Supply sum = 0; |
322 Supply sum = 0; |
397 for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n]; |
323 for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n]; |
398 valid_supply = sum == 0; |
324 valid_supply = sum == 0; |
410 /// \param _t The target node. |
336 /// \param _t The target node. |
411 /// \param _flow_value The required amount of flow from node \c _s |
337 /// \param _flow_value The required amount of flow from node \c _s |
412 /// to node \c _t (i.e. the supply of \c _s and the demand of |
338 /// to node \c _t (i.e. the supply of \c _s and the demand of |
413 /// \c _t). |
339 /// \c _t). |
414 CapacityScaling( const Graph &_graph, |
340 CapacityScaling( const Graph &_graph, |
415 const LowerMap &_lower, |
341 const LowerMap &_lower, |
416 const CapacityMap &_capacity, |
342 const CapacityMap &_capacity, |
417 const CostMap &_cost, |
343 const CostMap &_cost, |
418 Node _s, Node _t, |
344 Node _s, Node _t, |
419 Supply _flow_value ) : |
345 Supply _flow_value ) : |
420 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost), |
346 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost), |
421 supply(_graph), flow(_graph, 0), potential(_graph, 0), |
347 supply(_graph), flow(_graph, 0), potential(_graph, 0), |
422 res_graph(_graph, capacity, flow), |
348 res_cap(_graph), excess(_graph), pred(_graph), |
423 red_cost(res_graph, cost, potential), imbalance(_graph), |
349 dijkstra(graph, flow, res_cap, cost, excess, potential) |
424 #ifdef WITH_SCALING |
|
425 delta(0), delta_filter(res_graph, delta), |
|
426 dres_graph(res_graph, delta_filter), |
|
427 dijkstra(dres_graph, red_cost), pred(dres_graph), |
|
428 delta_deficit(imbalance, delta) |
|
429 #else //WITHOUT_CAPACITY_SCALING |
|
430 dijkstra(res_graph, red_cost), pred(res_graph), |
|
431 has_deficit(imbalance) |
|
432 #endif |
|
433 { |
350 { |
434 // Removing nonzero lower bounds |
351 // Removing nonzero lower bounds |
435 capacity = subMap(_capacity, _lower); |
352 capacity = subMap(_capacity, _lower); |
|
353 res_cap = capacity; |
436 for (NodeIt n(graph); n != INVALID; ++n) { |
354 for (NodeIt n(graph); n != INVALID; ++n) { |
437 Supply s = 0; |
355 Supply s = 0; |
438 if (n == _s) s = _flow_value; |
356 if (n == _s) s = _flow_value; |
439 if (n == _t) s = -_flow_value; |
357 if (n == _t) s = -_flow_value; |
440 for (InEdgeIt e(graph, n); e != INVALID; ++e) |
358 for (InEdgeIt e(graph, n); e != INVALID; ++e) |
441 s += _lower[e]; |
359 s += _lower[e]; |
442 for (OutEdgeIt e(graph, n); e != INVALID; ++e) |
360 for (OutEdgeIt e(graph, n); e != INVALID; ++e) |
443 s -= _lower[e]; |
361 s -= _lower[e]; |
444 supply[n] = s; |
362 supply[n] = s; |
445 } |
363 } |
446 valid_supply = true; |
364 valid_supply = true; |
447 } |
365 } |
448 |
366 |
449 /// \brief Simple constructor of the class (without lower bounds). |
367 /// \brief Simple constructor of the class (without lower bounds). |
457 /// \param _t The target node. |
375 /// \param _t The target node. |
458 /// \param _flow_value The required amount of flow from node \c _s |
376 /// \param _flow_value The required amount of flow from node \c _s |
459 /// to node \c _t (i.e. the supply of \c _s and the demand of |
377 /// to node \c _t (i.e. the supply of \c _s and the demand of |
460 /// \c _t). |
378 /// \c _t). |
461 CapacityScaling( const Graph &_graph, |
379 CapacityScaling( const Graph &_graph, |
462 const CapacityMap &_capacity, |
380 const CapacityMap &_capacity, |
463 const CostMap &_cost, |
381 const CostMap &_cost, |
464 Node _s, Node _t, |
382 Node _s, Node _t, |
465 Supply _flow_value ) : |
383 Supply _flow_value ) : |
466 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost), |
384 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost), |
467 supply(_graph, 0), flow(_graph, 0), potential(_graph, 0), |
385 supply(_graph, 0), flow(_graph, 0), potential(_graph, 0), |
468 res_graph(_graph, capacity, flow), |
386 res_cap(_capacity), excess(_graph), pred(_graph), |
469 red_cost(res_graph, cost, potential), imbalance(_graph), |
387 dijkstra(graph, flow, res_cap, cost, excess, potential) |
470 #ifdef WITH_SCALING |
|
471 delta(0), delta_filter(res_graph, delta), |
|
472 dres_graph(res_graph, delta_filter), |
|
473 dijkstra(dres_graph, red_cost), pred(dres_graph), |
|
474 delta_deficit(imbalance, delta) |
|
475 #else //WITHOUT_CAPACITY_SCALING |
|
476 dijkstra(res_graph, red_cost), pred(res_graph), |
|
477 has_deficit(imbalance) |
|
478 #endif |
|
479 { |
388 { |
480 supply[_s] = _flow_value; |
389 supply[_s] = _flow_value; |
481 supply[_t] = -_flow_value; |
390 supply[_t] = -_flow_value; |
482 valid_supply = true; |
391 valid_supply = true; |
483 } |
392 } |
509 /// |
418 /// |
510 /// \pre \ref run() must be called before using this function. |
419 /// \pre \ref run() must be called before using this function. |
511 Cost totalCost() const { |
420 Cost totalCost() const { |
512 Cost c = 0; |
421 Cost c = 0; |
513 for (EdgeIt e(graph); e != INVALID; ++e) |
422 for (EdgeIt e(graph); e != INVALID; ++e) |
514 c += flow[e] * cost[e]; |
423 c += flow[e] * cost[e]; |
515 return c; |
424 return c; |
516 } |
425 } |
517 |
426 |
518 /// \brief Runs the algorithm. |
427 /// \brief Runs the algorithm. |
519 /// |
428 /// |
520 /// Runs the algorithm. |
429 /// Runs the algorithm. |
521 /// |
430 /// |
|
431 /// \param scaling_mode The scaling mode. In case of WITH_SCALING |
|
432 /// capacity scaling is enabled in the algorithm (this is the |
|
433 /// default value) otherwise it is disabled. |
|
434 /// If the maximum edge capacity and/or the amount of total supply |
|
435 /// is small, the algorithm could be faster without scaling. |
|
436 /// |
522 /// \return \c true if a feasible flow can be found. |
437 /// \return \c true if a feasible flow can be found. |
523 bool run() { |
438 bool run(int scaling_mode = WITH_SCALING) { |
524 return init() && start(); |
439 return init(scaling_mode) && start(); |
525 } |
440 } |
526 |
441 |
527 protected: |
442 protected: |
528 |
443 |
529 /// \brief Initializes the algorithm. |
444 /// \brief Initializes the algorithm. |
530 bool init() { |
445 bool init(int scaling_mode) { |
531 if (!valid_supply) return false; |
446 if (!valid_supply) return false; |
532 imbalance = supply; |
447 excess = supply; |
533 |
448 |
534 // Initalizing Dijkstra class |
|
535 updater.potentialMap(potential); |
|
536 dijkstra.distMap(updater).predMap(pred); |
|
537 |
|
538 #ifdef WITH_SCALING |
|
539 // Initilaizing delta value |
449 // Initilaizing delta value |
540 Supply max_sup = 0, max_dem = 0; |
450 if (scaling_mode == WITH_SCALING) { |
541 for (NodeIt n(graph); n != INVALID; ++n) { |
451 // With scaling |
542 if (supply[n] > max_sup) max_sup = supply[n]; |
452 Supply max_sup = 0, max_dem = 0; |
543 if (supply[n] < -max_dem) max_dem = -supply[n]; |
453 for (NodeIt n(graph); n != INVALID; ++n) { |
544 } |
454 if ( supply[n] > max_sup) max_sup = supply[n]; |
545 if (max_dem < max_sup) max_sup = max_dem; |
455 if (-supply[n] > max_dem) max_dem = -supply[n]; |
546 for ( delta = 1; SCALING_FACTOR * delta < max_sup; |
456 } |
547 delta *= SCALING_FACTOR ) ; |
457 if (max_dem < max_sup) max_sup = max_dem; |
548 #endif |
458 phase_num = 0; |
|
459 for (delta = 1; 2 * delta <= max_sup; delta *= 2) |
|
460 ++phase_num; |
|
461 } else { |
|
462 // Without scaling |
|
463 delta = 1; |
|
464 } |
549 return true; |
465 return true; |
550 } |
466 } |
551 |
467 |
552 #ifdef WITH_SCALING |
468 /// \brief Executes the algorithm. |
|
469 bool start() { |
|
470 if (delta > 1) |
|
471 return startWithScaling(); |
|
472 else |
|
473 return startWithoutScaling(); |
|
474 } |
|
475 |
553 /// \brief Executes the capacity scaling version of the successive |
476 /// \brief Executes the capacity scaling version of the successive |
554 /// shortest path algorithm. |
477 /// shortest path algorithm. |
555 bool start() { |
478 bool startWithScaling() { |
556 typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt; |
|
557 typedef typename DeltaResGraph::Edge DeltaResEdge; |
|
558 #ifdef _DEBUG_ITER_ |
|
559 int dijk_num = 0; |
|
560 #endif |
|
561 |
|
562 // Processing capacity scaling phases |
479 // Processing capacity scaling phases |
563 ResNode s, t; |
480 Node s, t; |
564 for ( ; delta >= 1; delta = delta < SCALING_FACTOR && delta > 1 ? |
481 int phase_cnt = 0; |
565 1 : delta / SCALING_FACTOR ) |
482 int factor = 4; |
566 { |
483 while (true) { |
567 // Saturating edges not satisfying the optimality condition |
484 // Saturating all edges not satisfying the optimality condition |
568 Capacity r; |
485 for (EdgeIt e(graph); e != INVALID; ++e) { |
569 for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) { |
486 Node u = graph.source(e), v = graph.target(e); |
570 if (red_cost[e] < 0) { |
487 Cost c = cost[e] - potential[u] + potential[v]; |
571 res_graph.augment(e, r = res_graph.rescap(e)); |
488 if (c < 0 && res_cap[e] >= delta) { |
572 imbalance[dres_graph.source(e)] -= r; |
489 excess[u] -= res_cap[e]; |
573 imbalance[dres_graph.target(e)] += r; |
490 excess[v] += res_cap[e]; |
574 } |
491 flow[e] = capacity[e]; |
575 } |
492 res_cap[e] = 0; |
576 |
493 } |
577 // Finding excess nodes |
494 else if (c > 0 && flow[e] >= delta) { |
578 excess_nodes.clear(); |
495 excess[u] += flow[e]; |
579 deficit_nodes.clear(); |
496 excess[v] -= flow[e]; |
580 for (ResNodeIt n(res_graph); n != INVALID; ++n) { |
497 flow[e] = 0; |
581 if (imbalance[n] >= delta) excess_nodes.push_back(n); |
498 res_cap[e] = capacity[e]; |
582 if (imbalance[n] <= -delta) deficit_nodes.push_back(n); |
499 } |
583 } |
500 } |
584 next_node = 0; |
501 |
585 |
502 // Finding excess nodes and deficit nodes |
586 // Finding shortest paths |
503 excess_nodes.clear(); |
587 while (next_node < excess_nodes.size()) { |
504 deficit_nodes.clear(); |
588 // Checking deficit nodes |
505 for (NodeIt n(graph); n != INVALID; ++n) { |
589 if (delta > 1) { |
506 if (excess[n] >= delta) excess_nodes.push_back(n); |
590 bool delta_def = false; |
507 if (excess[n] <= -delta) deficit_nodes.push_back(n); |
591 for (int i = 0; i < deficit_nodes.size(); ++i) { |
508 } |
592 if (imbalance[deficit_nodes[i]] <= -delta) { |
509 next_node = 0; |
593 delta_def = true; |
510 |
594 break; |
511 // Finding augmenting shortest paths |
595 } |
512 while (next_node < excess_nodes.size()) { |
596 } |
513 // Checking deficit nodes |
597 if (!delta_def) break; |
514 if (delta > 1) { |
598 } |
515 bool delta_deficit = false; |
599 |
516 for (int i = 0; i < deficit_nodes.size(); ++i) { |
600 // Running Dijkstra |
517 if (excess[deficit_nodes[i]] <= -delta) { |
601 s = excess_nodes[next_node]; |
518 delta_deficit = true; |
602 updater.init(); |
519 break; |
603 dijkstra.init(); |
520 } |
604 dijkstra.addSource(s); |
521 } |
605 #ifdef _DEBUG_ITER_ |
522 if (!delta_deficit) break; |
606 ++dijk_num; |
523 } |
607 #endif |
524 |
608 if ((t = dijkstra.start(delta_deficit)) == INVALID) { |
525 // Running Dijkstra |
609 if (delta > 1) { |
526 s = excess_nodes[next_node]; |
610 ++next_node; |
527 if ((t = dijkstra.run(s, delta)) == INVALID) { |
611 continue; |
528 if (delta > 1) { |
612 } |
529 ++next_node; |
613 return false; |
530 continue; |
614 } |
531 } |
615 |
532 return false; |
616 // Updating node potentials |
533 } |
617 updater.update(); |
534 |
618 |
535 // Augmenting along a shortest path from s to t. |
619 // Augment along a shortest path from s to t |
536 Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t]; |
620 Capacity d = imbalance[s] < -imbalance[t] ? |
537 Node u = t; |
621 imbalance[s] : -imbalance[t]; |
538 Edge e; |
622 ResNode u = t; |
539 if (d > delta) { |
623 ResEdge e; |
540 while ((e = pred[u]) != INVALID) { |
624 if (d > delta) { |
541 Capacity rc; |
625 while ((e = pred[u]) != INVALID) { |
542 if (u == graph.target(e)) { |
626 if (res_graph.rescap(e) < d) d = res_graph.rescap(e); |
543 rc = res_cap[e]; |
627 u = dres_graph.source(e); |
544 u = graph.source(e); |
628 } |
545 } else { |
629 } |
546 rc = flow[e]; |
630 u = t; |
547 u = graph.target(e); |
631 while ((e = pred[u]) != INVALID) { |
548 } |
632 res_graph.augment(e, d); |
549 if (rc < d) d = rc; |
633 u = dres_graph.source(e); |
550 } |
634 } |
551 } |
635 imbalance[s] -= d; |
552 u = t; |
636 imbalance[t] += d; |
553 while ((e = pred[u]) != INVALID) { |
637 if (imbalance[s] < delta) ++next_node; |
554 if (u == graph.target(e)) { |
638 } |
555 flow[e] += d; |
639 } |
556 res_cap[e] -= d; |
640 #ifdef _DEBUG_ITER_ |
557 u = graph.source(e); |
641 std::cout << "Capacity Scaling algorithm finished with running Dijkstra algorithm " |
558 } else { |
642 << dijk_num << " times." << std::endl; |
559 flow[e] -= d; |
643 #endif |
560 res_cap[e] += d; |
|
561 u = graph.target(e); |
|
562 } |
|
563 } |
|
564 excess[s] -= d; |
|
565 excess[t] += d; |
|
566 |
|
567 if (excess[s] < delta) ++next_node; |
|
568 } |
|
569 |
|
570 if (delta == 1) break; |
|
571 if (++phase_cnt > phase_num / 4) factor = 2; |
|
572 delta = delta <= factor ? 1 : delta / factor; |
|
573 } |
644 |
574 |
645 // Handling nonzero lower bounds |
575 // Handling nonzero lower bounds |
646 if (lower) { |
576 if (lower) { |
647 for (EdgeIt e(graph); e != INVALID; ++e) |
577 for (EdgeIt e(graph); e != INVALID; ++e) |
648 flow[e] += (*lower)[e]; |
578 flow[e] += (*lower)[e]; |
649 } |
579 } |
650 return true; |
580 return true; |
651 } |
581 } |
652 |
582 |
653 #else //WITHOUT_CAPACITY_SCALING |
|
654 /// \brief Executes the successive shortest path algorithm without |
583 /// \brief Executes the successive shortest path algorithm without |
655 /// capacity scaling. |
584 /// capacity scaling. |
656 bool start() { |
585 bool startWithoutScaling() { |
657 // Finding excess nodes |
586 // Finding excess nodes |
658 for (ResNodeIt n(res_graph); n != INVALID; ++n) { |
587 for (NodeIt n(graph); n != INVALID; ++n) { |
659 if (imbalance[n] > 0) excess_nodes.push_back(n); |
588 if (excess[n] > 0) excess_nodes.push_back(n); |
660 } |
589 } |
661 if (excess_nodes.size() == 0) return true; |
590 if (excess_nodes.size() == 0) return true; |
662 next_node = 0; |
591 next_node = 0; |
663 |
592 |
664 // Finding shortest paths |
593 // Finding shortest paths |
665 ResNode s, t; |
594 Node s, t; |
666 while ( imbalance[excess_nodes[next_node]] > 0 || |
595 while ( excess[excess_nodes[next_node]] > 0 || |
667 ++next_node < excess_nodes.size() ) |
596 ++next_node < excess_nodes.size() ) |
668 { |
597 { |
669 // Running Dijkstra |
598 // Running Dijkstra |
670 s = excess_nodes[next_node]; |
599 s = excess_nodes[next_node]; |
671 updater.init(); |
600 if ((t = dijkstra.run(s, 1)) == INVALID) |
672 dijkstra.init(); |
601 return false; |
673 dijkstra.addSource(s); |
602 |
674 if ((t = dijkstra.start(has_deficit)) == INVALID) |
603 // Augmenting along a shortest path from s to t |
675 return false; |
604 Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t]; |
676 |
605 Node u = t; |
677 // Updating node potentials |
606 Edge e; |
678 updater.update(); |
607 while ((e = pred[u]) != INVALID) { |
679 |
608 Capacity rc; |
680 // Augmenting along a shortest path from s to t |
609 if (u == graph.target(e)) { |
681 Capacity delta = imbalance[s] < -imbalance[t] ? |
610 rc = res_cap[e]; |
682 imbalance[s] : -imbalance[t]; |
611 u = graph.source(e); |
683 ResNode u = t; |
612 } else { |
684 ResEdge e; |
613 rc = flow[e]; |
685 while ((e = pred[u]) != INVALID) { |
614 u = graph.target(e); |
686 if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e); |
615 } |
687 u = res_graph.source(e); |
616 if (rc < d) d = rc; |
688 } |
617 } |
689 u = t; |
618 u = t; |
690 while ((e = pred[u]) != INVALID) { |
619 while ((e = pred[u]) != INVALID) { |
691 res_graph.augment(e, delta); |
620 if (u == graph.target(e)) { |
692 u = res_graph.source(e); |
621 flow[e] += d; |
693 } |
622 res_cap[e] -= d; |
694 imbalance[s] -= delta; |
623 u = graph.source(e); |
695 imbalance[t] += delta; |
624 } else { |
|
625 flow[e] -= d; |
|
626 res_cap[e] += d; |
|
627 u = graph.target(e); |
|
628 } |
|
629 } |
|
630 excess[s] -= d; |
|
631 excess[t] += d; |
696 } |
632 } |
697 |
633 |
698 // Handling nonzero lower bounds |
634 // Handling nonzero lower bounds |
699 if (lower) { |
635 if (lower) { |
700 for (EdgeIt e(graph); e != INVALID; ++e) |
636 for (EdgeIt e(graph); e != INVALID; ++e) |
701 flow[e] += (*lower)[e]; |
637 flow[e] += (*lower)[e]; |
702 } |
638 } |
703 return true; |
639 return true; |
704 } |
640 } |
705 #endif |
|
706 |
641 |
707 }; //class CapacityScaling |
642 }; //class CapacityScaling |
708 |
643 |
709 ///@} |
644 ///@} |
710 |
645 |