A sample file.
2 #ifndef LEMON_MINCOSTFLOW_H
3 #define LEMON_MINCOSTFLOW_H
7 ///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network
9 #include <lemon/dijkstra.h>
10 #include <lemon/graph_wrapper.h>
11 #include <lemon/maps.h>
15 #include <lemon/for_each_macros.h>
16 #include <lemon/unionfind.h>
17 #include <lemon/bin_heap.h>
25 ///\brief Implementation of an algorithm for solving the minimum cost general
26 /// flow problem in an uncapacitated network
29 /// The class \ref lemon::MinCostFlow "MinCostFlow" implements
30 /// an algorithm for solving the following general minimum cost flow problem>
34 /// \warning It is assumed here that the problem has a feasible solution
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.
42 ///\author Attila Bernath
43 template <typename Graph, typename CostMap, typename SupplyDemandMap>
46 typedef typename CostMap::Value Cost;
49 typedef typename SupplyDemandMap::Value SupplyDemand;
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;
58 // typedef ConstMap<Edge,int> ConstMap;
60 typedef ResGraphWrapper<const Graph,int,ConstEdgeMap,FlowMap> ResGraph;
61 typedef typename ResGraph::Edge ResGraphEdge;
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;
71 typedef typename CostMap::Key Key;
72 typedef typename CostMap::Value Value;
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)]);
78 return -ol[e]-(pot[res_graph.target(e)]-pot[res_graph.source(e)]);
81 ModCostMap(const ResGraph& _res_graph,
82 const CostMap &o, const NodeMap &p) :
83 res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
92 const SupplyDemandMap& supply_demand;//supply or demand of nodes
99 //To store the potential (dual variables)
100 typedef typename Graph::template NodeMap<Cost> PotentialMap;
101 PotentialMap potential;
110 MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand):
113 supply_demand(_supply_demand),
118 ///Runs the algorithm.
120 ///Runs the algorithm.
122 ///\todo May be it does make sense to be able to start with a nonzero
123 /// feasible primal-dual solution pair as well.
126 //To store excess-deficit values
127 SupplyDemandMap excess_deficit(graph);
129 //Resetting variables from previous runs
133 typedef typename Graph::template NodeMap<int> HeapMap;
134 typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
135 std::greater<SupplyDemand> > HeapType;
137 //A heap for the excess nodes
138 HeapMap excess_nodes_map(graph,-1);
139 HeapType excess_nodes(excess_nodes_map);
141 //A heap for the deficit nodes
142 HeapMap deficit_nodes_map(graph,-1);
143 HeapType deficit_nodes(deficit_nodes_map);
145 //A container to store nonabundant arcs
146 std::list<Edge> nonabundant_arcs;
149 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
151 nonabundant_arcs.push_back(e);
154 //Initial value for delta
155 SupplyDemand delta = 0;
157 typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
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);
165 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
166 excess_deficit.set(n,supply_demand[n]);
168 if (excess_deficit[n] > 0){
169 excess_nodes.push(n,excess_deficit[n]);
172 if (excess_deficit[n] < 0){
173 deficit_nodes.push(n, - excess_deficit[n]);
175 //Finding out starting value of delta
176 if (delta < abs(excess_deficit[n])){
177 delta = abs(excess_deficit[n]);
179 //Initialize the copy of the Dijkstra potential to zero
181 //Every single point is an abundant component initially
182 abundant_components.insert(n);
185 //It'll be allright as an initial value, though this value
186 //can be the maximum deficit here
187 SupplyDemand max_excess = delta;
189 ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph
190 ConstEdgeMap const_inf_map(MAXINT);
192 //We need a residual graph which is uncapacitated
193 ResGraph res_graph(graph, const_inf_map, flow);
195 //An EdgeMap to tell which arcs are abundant
196 typename Graph::template EdgeMap<bool> abundant_arcs(graph);
198 //Let's construct the sugraph consisting only of the abundant edges
199 typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
201 ConstNodeMap const_true_map(true);
202 typedef SubGraphWrapper< const Graph, ConstNodeMap,
203 typename Graph::template EdgeMap<bool> >
205 AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
207 //Let's construct the residual graph for the abundant graph
208 typedef ResGraphWrapper<const AbundantGraph,int,ConstEdgeMap,FlowMap>
210 //Again uncapacitated
211 ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow);
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;
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>,
231 =typename Graph::template NodeMap<typename Graph::Edge>,
232 typename DistMap=typename Graph::template NodeMap<int> >
233 class Bfs : public BfsIterator<Graph, ReachedMap> {
237 ModCostMap mod_cost(res_graph, cost, potential);
239 Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);
241 //We will use the number of the nodes of the graph often
242 int number_of_nodes = graph.nodeNum();
244 while (max_excess > 0){
246 //Reset delta if still too big
247 if (8*number_of_nodes*max_excess <= delta){
253 * Beginning of the delta scaling phase
257 SupplyDemand buf=8*number_of_nodes*delta;
258 typename std::list<Edge>::iterator i = nonabundant_arcs.begin();
259 while ( i != nonabundant_arcs.end() ){
261 Node a = abundant_components.find(res_graph.target(*i));
262 Node b = abundant_components.find(res_graph.source(*i));
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;
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
282 //Augmenting on the found path
287 n = res_graph.source(e);
288 res_graph.augment(e,qty_to_augment);
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);
301 //Or negative but not turned into positive
302 deficit_nodes.set(root,
303 deficit_nodes[root] - qty_to_augment);
306 //Update the excess_deficit map
307 excess_deficit[non_root] -= qty_to_augment;
308 excess_deficit[root] += qty_to_augment;
314 //Marci and Zsolt says I shouldn't do such things
315 nonabundant_arcs.erase(i++);
316 abundant_arcs[*i] = true;
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];
332 while(max_excess > (number_of_nodes-1)*delta/number_of_nodes){
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
347 //We have to change the potential
348 FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
349 potential[n] += dijkstra.distMap()[n];
353 //Augmenting on the sortest path
357 e = dijkstra.pred(n);
358 n = dijkstra.predNode(n);
359 res_graph.augment(e,delta);
361 //Let's update the total cost
362 if (res_graph.forward(e))
363 total_cost += cost[e];
365 total_cost -= cost[e];
369 //Update the excess_deficit map
370 excess_deficit[s] -= delta;
371 excess_deficit[t] += delta;
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]);
382 excess_nodes.set(s, excess_nodes[s] - delta);
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]);
392 deficit_nodes.set(t, deficit_nodes[t] - delta);
394 //Dijkstra part ends here
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];
407 * End of the delta scaling phase
410 //Whatever this means
413 /*This is not necessary here
414 //Update the max_excess
416 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
417 if (max_excess < excess_deficit[n]){
418 max_excess = excess_deficit[n];
424 }//while(max_excess > 0)
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.
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;}
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;}
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
451 ///\todo Is this OK here?
452 bool checkComplementarySlackness(){
455 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
457 mod_pot = cost[e]-potential[graph.target(e)]+potential[graph.source(e)];
459 // std::cout << fl_e << std::endl;
460 if (mod_pot > 0 && fl_e != 0)
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;
478 bool checkFeasibility(){
479 SupplyDemandMap supdem(graph);
480 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
486 supdem[graph.source(e)] += flow[e];
487 supdem[graph.target(e)] -= flow[e];
489 //write_property_vector(supdem, "supdem");
490 //write_property_vector(supply_demand, "supply_demand");
492 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
494 if ( supdem[n] != supply_demand[n]){
495 //cout<<"Node id.: "<<graph.id(n)<<" : "<<supdem[n]<<", should be: "<<supply_demand[n]<<endl;
503 bool checkOptimality(){
504 return checkFeasibility() && checkComplementarySlackness();
507 }; //class MinCostFlow
513 #endif //LEMON_MINCOSTFLOW_H