1 // -*- c++ -*- |
|
2 #ifndef LEMON_MINCOSTFLOW_H |
|
3 #define LEMON_MINCOSTFLOW_H |
|
4 |
|
5 ///\ingroup galgs |
|
6 ///\file |
|
7 ///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network |
|
8 |
|
9 #include <lemon/dijkstra.h> |
|
10 #include <lemon/graph_wrapper.h> |
|
11 #include <lemon/maps.h> |
|
12 #include <vector> |
|
13 #include <list> |
|
14 #include <values.h> |
|
15 #include <lemon/for_each_macros.h> |
|
16 #include <lemon/unionfind.h> |
|
17 #include <lemon/bin_heap.h> |
|
18 #include <bfs_dfs.h> |
|
19 |
|
20 namespace lemon { |
|
21 |
|
22 /// \addtogroup galgs |
|
23 /// @{ |
|
24 |
|
25 ///\brief Implementation of an algorithm for solving the minimum cost general |
|
26 /// flow problem in an uncapacitated network |
|
27 /// |
|
28 /// |
|
29 /// The class \ref lemon::MinCostFlow "MinCostFlow" implements |
|
30 /// an algorithm for solving the following general minimum cost flow problem> |
|
31 /// |
|
32 /// |
|
33 /// |
|
34 /// \warning It is assumed here that the problem has a feasible solution |
|
35 /// |
|
36 /// The range of the cost (weight) function is nonnegative reals but |
|
37 /// the range of capacity function is the set of nonnegative integers. |
|
38 /// It is not a polinomial time algorithm for counting the minimum cost |
|
39 /// maximal flow, since it counts the minimum cost flow for every value 0..M |
|
40 /// where \c M is the value of the maximal flow. |
|
41 /// |
|
42 ///\author Attila Bernath |
|
43 template <typename Graph, typename CostMap, typename SupplyDemandMap> |
|
44 class MinCostFlow { |
|
45 |
|
46 typedef typename CostMap::Value Cost; |
|
47 |
|
48 |
|
49 typedef typename SupplyDemandMap::Value SupplyDemand; |
|
50 |
|
51 typedef typename Graph::Node Node; |
|
52 typedef typename Graph::NodeIt NodeIt; |
|
53 typedef typename Graph::Edge Edge; |
|
54 typedef typename Graph::OutEdgeIt OutEdgeIt; |
|
55 typedef typename Graph::template EdgeMap<SupplyDemand> FlowMap; |
|
56 typedef ConstMap<Edge,SupplyDemand> ConstEdgeMap; |
|
57 |
|
58 // typedef ConstMap<Edge,int> ConstMap; |
|
59 |
|
60 typedef ResGraphWrapper<const Graph,int,ConstEdgeMap,FlowMap> ResGraph; |
|
61 typedef typename ResGraph::Edge ResGraphEdge; |
|
62 |
|
63 class ModCostMap { |
|
64 //typedef typename ResGraph::template NodeMap<Cost> NodeMap; |
|
65 typedef typename Graph::template NodeMap<Cost> NodeMap; |
|
66 const ResGraph& res_graph; |
|
67 // const EdgeIntMap& rev; |
|
68 const CostMap &ol; |
|
69 const NodeMap &pot; |
|
70 public : |
|
71 typedef typename CostMap::Key Key; |
|
72 typedef typename CostMap::Value Value; |
|
73 |
|
74 Value operator[](typename ResGraph::Edge e) const { |
|
75 if (res_graph.forward(e)) |
|
76 return ol[e]-(pot[res_graph.target(e)]-pot[res_graph.source(e)]); |
|
77 else |
|
78 return -ol[e]-(pot[res_graph.target(e)]-pot[res_graph.source(e)]); |
|
79 } |
|
80 |
|
81 ModCostMap(const ResGraph& _res_graph, |
|
82 const CostMap &o, const NodeMap &p) : |
|
83 res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){}; |
|
84 };//ModCostMap |
|
85 |
|
86 |
|
87 protected: |
|
88 |
|
89 //Input |
|
90 const Graph& graph; |
|
91 const CostMap& cost; |
|
92 const SupplyDemandMap& supply_demand;//supply or demand of nodes |
|
93 |
|
94 |
|
95 //auxiliary variables |
|
96 |
|
97 //To store the flow |
|
98 FlowMap flow; |
|
99 //To store the potential (dual variables) |
|
100 typedef typename Graph::template NodeMap<Cost> PotentialMap; |
|
101 PotentialMap potential; |
|
102 |
|
103 |
|
104 Cost total_cost; |
|
105 |
|
106 |
|
107 public : |
|
108 |
|
109 |
|
110 MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand): |
|
111 graph(_graph), |
|
112 cost(_cost), |
|
113 supply_demand(_supply_demand), |
|
114 flow(_graph), |
|
115 potential(_graph){ } |
|
116 |
|
117 |
|
118 ///Runs the algorithm. |
|
119 |
|
120 ///Runs the algorithm. |
|
121 |
|
122 ///\todo May be it does make sense to be able to start with a nonzero |
|
123 /// feasible primal-dual solution pair as well. |
|
124 void run() { |
|
125 |
|
126 //To store excess-deficit values |
|
127 SupplyDemandMap excess_deficit(graph); |
|
128 |
|
129 //Resetting variables from previous runs |
|
130 //total_cost = 0; |
|
131 |
|
132 |
|
133 typedef typename Graph::template NodeMap<int> HeapMap; |
|
134 typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap<int>, |
|
135 std::greater<SupplyDemand> > HeapType; |
|
136 |
|
137 //A heap for the excess nodes |
|
138 HeapMap excess_nodes_map(graph,-1); |
|
139 HeapType excess_nodes(excess_nodes_map); |
|
140 |
|
141 //A heap for the deficit nodes |
|
142 HeapMap deficit_nodes_map(graph,-1); |
|
143 HeapType deficit_nodes(deficit_nodes_map); |
|
144 |
|
145 //A container to store nonabundant arcs |
|
146 std::list<Edge> nonabundant_arcs; |
|
147 |
|
148 |
|
149 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ |
|
150 flow.set(e,0); |
|
151 nonabundant_arcs.push_back(e); |
|
152 } |
|
153 |
|
154 //Initial value for delta |
|
155 SupplyDemand delta = 0; |
|
156 |
|
157 typedef UnionFindEnum<Node, Graph::template NodeMap> UFE; |
|
158 |
|
159 //A union-find structure to store the abundant components |
|
160 typename UFE::MapType abund_comp_map(graph); |
|
161 UFE abundant_components(abund_comp_map); |
|
162 |
|
163 |
|
164 |
|
165 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){ |
|
166 excess_deficit.set(n,supply_demand[n]); |
|
167 //A supply node |
|
168 if (excess_deficit[n] > 0){ |
|
169 excess_nodes.push(n,excess_deficit[n]); |
|
170 } |
|
171 //A demand node |
|
172 if (excess_deficit[n] < 0){ |
|
173 deficit_nodes.push(n, - excess_deficit[n]); |
|
174 } |
|
175 //Finding out starting value of delta |
|
176 if (delta < abs(excess_deficit[n])){ |
|
177 delta = abs(excess_deficit[n]); |
|
178 } |
|
179 //Initialize the copy of the Dijkstra potential to zero |
|
180 potential.set(n,0); |
|
181 //Every single point is an abundant component initially |
|
182 abundant_components.insert(n); |
|
183 } |
|
184 |
|
185 //It'll be allright as an initial value, though this value |
|
186 //can be the maximum deficit here |
|
187 SupplyDemand max_excess = delta; |
|
188 |
|
189 ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph |
|
190 ConstEdgeMap const_inf_map(MAXINT); |
|
191 |
|
192 //We need a residual graph which is uncapacitated |
|
193 ResGraph res_graph(graph, const_inf_map, flow); |
|
194 |
|
195 //An EdgeMap to tell which arcs are abundant |
|
196 typename Graph::template EdgeMap<bool> abundant_arcs(graph); |
|
197 |
|
198 //Let's construct the sugraph consisting only of the abundant edges |
|
199 typedef ConstMap< typename Graph::Node, bool > ConstNodeMap; |
|
200 |
|
201 ConstNodeMap const_true_map(true); |
|
202 typedef SubGraphWrapper< const Graph, ConstNodeMap, |
|
203 typename Graph::template EdgeMap<bool> > |
|
204 AbundantGraph; |
|
205 AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs ); |
|
206 |
|
207 //Let's construct the residual graph for the abundant graph |
|
208 typedef ResGraphWrapper<const AbundantGraph,int,ConstEdgeMap,FlowMap> |
|
209 ResAbGraph; |
|
210 //Again uncapacitated |
|
211 ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow); |
|
212 |
|
213 //We need things for the bfs |
|
214 typename ResAbGraph::template NodeMap<bool> bfs_reached(res_ab_graph); |
|
215 typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge> |
|
216 bfs_pred(res_ab_graph); |
|
217 NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy; |
|
218 //Teszt celbol: |
|
219 //BfsIterator<ResAbGraph, typename ResAbGraph::template NodeMap<bool> > |
|
220 //izebize(res_ab_graph, bfs_reached); |
|
221 //We want to run bfs-es (more) on this graph 'res_ab_graph' |
|
222 Bfs < const ResAbGraph , |
|
223 typename ResAbGraph::template NodeMap<bool>, |
|
224 typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>, |
|
225 NullMap<typename ResAbGraph::Node, int> > |
|
226 bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy); |
|
227 /*This is what Marci wants for a bfs |
|
228 template <typename Graph, |
|
229 typename ReachedMap=typename Graph::template NodeMap<bool>, |
|
230 typename PredMap |
|
231 =typename Graph::template NodeMap<typename Graph::Edge>, |
|
232 typename DistMap=typename Graph::template NodeMap<int> > |
|
233 class Bfs : public BfsIterator<Graph, ReachedMap> { |
|
234 |
|
235 */ |
|
236 |
|
237 ModCostMap mod_cost(res_graph, cost, potential); |
|
238 |
|
239 Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost); |
|
240 |
|
241 //We will use the number of the nodes of the graph often |
|
242 int number_of_nodes = graph.nodeNum(); |
|
243 |
|
244 while (max_excess > 0){ |
|
245 |
|
246 //Reset delta if still too big |
|
247 if (8*number_of_nodes*max_excess <= delta){ |
|
248 delta = max_excess; |
|
249 |
|
250 } |
|
251 |
|
252 /* |
|
253 * Beginning of the delta scaling phase |
|
254 */ |
|
255 //Merge and stuff |
|
256 { |
|
257 SupplyDemand buf=8*number_of_nodes*delta; |
|
258 typename std::list<Edge>::iterator i = nonabundant_arcs.begin(); |
|
259 while ( i != nonabundant_arcs.end() ){ |
|
260 if (flow[*i]>=buf){ |
|
261 Node a = abundant_components.find(res_graph.target(*i)); |
|
262 Node b = abundant_components.find(res_graph.source(*i)); |
|
263 //Merge |
|
264 if (a != b){ |
|
265 abundant_components.join(a,b); |
|
266 //We want to push the smaller |
|
267 //Which has greater absolut value excess/deficit |
|
268 Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b; |
|
269 //Which is the other |
|
270 Node non_root = ( a == root ) ? b : a ; |
|
271 abundant_components.makeRep(root); |
|
272 SupplyDemand qty_to_augment = abs(excess_deficit[non_root]); |
|
273 //Push the positive value |
|
274 if (excess_deficit[non_root] < 0) |
|
275 swap(root, non_root); |
|
276 //If the non_root node has excess/deficit at all |
|
277 if (qty_to_augment>0){ |
|
278 //Find path and augment |
|
279 bfs.run(typename AbundantGraph::Node(non_root)); |
|
280 //root should be reached |
|
281 |
|
282 //Augmenting on the found path |
|
283 Node n=root; |
|
284 ResGraphEdge e; |
|
285 while (n!=non_root){ |
|
286 e = bfs_pred[n]; |
|
287 n = res_graph.source(e); |
|
288 res_graph.augment(e,qty_to_augment); |
|
289 } |
|
290 |
|
291 //We know that non_root had positive excess |
|
292 excess_nodes.set(non_root, |
|
293 excess_nodes[non_root] - qty_to_augment); |
|
294 //But what about root node |
|
295 //It might have been positive and so became larger |
|
296 if (excess_deficit[root]>0){ |
|
297 excess_nodes.set(root, |
|
298 excess_nodes[root] + qty_to_augment); |
|
299 } |
|
300 else{ |
|
301 //Or negative but not turned into positive |
|
302 deficit_nodes.set(root, |
|
303 deficit_nodes[root] - qty_to_augment); |
|
304 } |
|
305 |
|
306 //Update the excess_deficit map |
|
307 excess_deficit[non_root] -= qty_to_augment; |
|
308 excess_deficit[root] += qty_to_augment; |
|
309 |
|
310 |
|
311 } |
|
312 } |
|
313 //What happens to i? |
|
314 //Marci and Zsolt says I shouldn't do such things |
|
315 nonabundant_arcs.erase(i++); |
|
316 abundant_arcs[*i] = true; |
|
317 } |
|
318 else |
|
319 ++i; |
|
320 } |
|
321 } |
|
322 |
|
323 |
|
324 Node s = excess_nodes.top(); |
|
325 max_excess = excess_nodes[s]; |
|
326 Node t = deficit_nodes.top(); |
|
327 if (max_excess < deficit_nodes[t]){ |
|
328 max_excess = deficit_nodes[t]; |
|
329 } |
|
330 |
|
331 |
|
332 while(max_excess > (number_of_nodes-1)*delta/number_of_nodes){ |
|
333 |
|
334 |
|
335 //s es t valasztasa |
|
336 |
|
337 //Dijkstra part |
|
338 dijkstra.run(s); |
|
339 |
|
340 /*We know from theory that t can be reached |
|
341 if (!dijkstra.reached(t)){ |
|
342 //There are no k paths from s to t |
|
343 break; |
|
344 }; |
|
345 */ |
|
346 |
|
347 //We have to change the potential |
|
348 FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){ |
|
349 potential[n] += dijkstra.distMap()[n]; |
|
350 } |
|
351 |
|
352 |
|
353 //Augmenting on the sortest path |
|
354 Node n=t; |
|
355 ResGraphEdge e; |
|
356 while (n!=s){ |
|
357 e = dijkstra.pred(n); |
|
358 n = dijkstra.predNode(n); |
|
359 res_graph.augment(e,delta); |
|
360 /* |
|
361 //Let's update the total cost |
|
362 if (res_graph.forward(e)) |
|
363 total_cost += cost[e]; |
|
364 else |
|
365 total_cost -= cost[e]; |
|
366 */ |
|
367 } |
|
368 |
|
369 //Update the excess_deficit map |
|
370 excess_deficit[s] -= delta; |
|
371 excess_deficit[t] += delta; |
|
372 |
|
373 |
|
374 //Update the excess_nodes heap |
|
375 if (delta > excess_nodes[s]){ |
|
376 if (delta > excess_nodes[s]) |
|
377 deficit_nodes.push(s,delta - excess_nodes[s]); |
|
378 excess_nodes.pop(); |
|
379 |
|
380 } |
|
381 else{ |
|
382 excess_nodes.set(s, excess_nodes[s] - delta); |
|
383 } |
|
384 //Update the deficit_nodes heap |
|
385 if (delta > deficit_nodes[t]){ |
|
386 if (delta > deficit_nodes[t]) |
|
387 excess_nodes.push(t,delta - deficit_nodes[t]); |
|
388 deficit_nodes.pop(); |
|
389 |
|
390 } |
|
391 else{ |
|
392 deficit_nodes.set(t, deficit_nodes[t] - delta); |
|
393 } |
|
394 //Dijkstra part ends here |
|
395 |
|
396 //Choose s and t again |
|
397 s = excess_nodes.top(); |
|
398 max_excess = excess_nodes[s]; |
|
399 t = deficit_nodes.top(); |
|
400 if (max_excess < deficit_nodes[t]){ |
|
401 max_excess = deficit_nodes[t]; |
|
402 } |
|
403 |
|
404 } |
|
405 |
|
406 /* |
|
407 * End of the delta scaling phase |
|
408 */ |
|
409 |
|
410 //Whatever this means |
|
411 delta = delta / 2; |
|
412 |
|
413 /*This is not necessary here |
|
414 //Update the max_excess |
|
415 max_excess = 0; |
|
416 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){ |
|
417 if (max_excess < excess_deficit[n]){ |
|
418 max_excess = excess_deficit[n]; |
|
419 } |
|
420 } |
|
421 */ |
|
422 |
|
423 |
|
424 }//while(max_excess > 0) |
|
425 |
|
426 |
|
427 //return i; |
|
428 } |
|
429 |
|
430 |
|
431 |
|
432 |
|
433 ///This function gives back the total cost of the found paths. |
|
434 ///Assumes that \c run() has been run and nothing changed since then. |
|
435 Cost totalCost(){ |
|
436 return total_cost; |
|
437 } |
|
438 |
|
439 ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must |
|
440 ///be called before using this function. |
|
441 const FlowMap &getFlow() const { return flow;} |
|
442 |
|
443 ///Returns a const reference to the NodeMap \c potential (the dual solution). |
|
444 /// \pre \ref run() must be called before using this function. |
|
445 const PotentialMap &getPotential() const { return potential;} |
|
446 |
|
447 ///This function checks, whether the given solution is optimal |
|
448 ///Running after a \c run() should return with true |
|
449 ///In this "state of the art" this only checks optimality, doesn't bother with feasibility |
|
450 /// |
|
451 ///\todo Is this OK here? |
|
452 bool checkComplementarySlackness(){ |
|
453 Cost mod_pot; |
|
454 Cost fl_e; |
|
455 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ |
|
456 //C^{\Pi}_{i,j} |
|
457 mod_pot = cost[e]-potential[graph.target(e)]+potential[graph.source(e)]; |
|
458 fl_e = flow[e]; |
|
459 // std::cout << fl_e << std::endl; |
|
460 if (mod_pot > 0 && fl_e != 0) |
|
461 return false; |
|
462 |
|
463 } |
|
464 return true; |
|
465 } |
|
466 |
|
467 /* |
|
468 //For testing purposes only |
|
469 //Lists the node_properties |
|
470 void write_property_vector(const SupplyDemandMap& a, |
|
471 char* prop_name="property"){ |
|
472 FOR_EACH_LOC(typename Graph::NodeIt, i, graph){ |
|
473 cout<<"Node id.: "<<graph.id(i)<<", "<<prop_name<<" value: "<<a[i]<<endl; |
|
474 } |
|
475 cout<<endl; |
|
476 } |
|
477 */ |
|
478 bool checkFeasibility(){ |
|
479 SupplyDemandMap supdem(graph); |
|
480 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ |
|
481 |
|
482 if ( flow[e] < 0){ |
|
483 |
|
484 return false; |
|
485 } |
|
486 supdem[graph.source(e)] += flow[e]; |
|
487 supdem[graph.target(e)] -= flow[e]; |
|
488 } |
|
489 //write_property_vector(supdem, "supdem"); |
|
490 //write_property_vector(supply_demand, "supply_demand"); |
|
491 |
|
492 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){ |
|
493 |
|
494 if ( supdem[n] != supply_demand[n]){ |
|
495 //cout<<"Node id.: "<<graph.id(n)<<" : "<<supdem[n]<<", should be: "<<supply_demand[n]<<endl; |
|
496 return false; |
|
497 } |
|
498 } |
|
499 |
|
500 return true; |
|
501 } |
|
502 |
|
503 bool checkOptimality(){ |
|
504 return checkFeasibility() && checkComplementarySlackness(); |
|
505 } |
|
506 |
|
507 }; //class MinCostFlow |
|
508 |
|
509 ///@} |
|
510 |
|
511 } //namespace lemon |
|
512 |
|
513 #endif //LEMON_MINCOSTFLOW_H |
|