Almost ready.
2 #ifndef HUGO_MINCOSTFLOW_H
3 #define HUGO_MINCOSTFLOW_H
7 ///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network
9 #include <hugo/dijkstra.h>
10 #include <hugo/graph_wrapper.h>
11 #include <hugo/maps.h>
14 #include <for_each_macros.h>
15 #include <hugo/union_find.h>
22 ///\brief Implementation of an algorithm for finding the minimum cost flow
23 /// of given value in an uncapacitated network
26 /// The class \ref hugo::MinCostFlow "MinCostFlow" implements
27 /// an algorithm for solving the following general minimum cost flow problem>
31 /// \warning It is assumed here that the problem has a feasible solution
33 /// The range of the length (weight) function is nonnegative reals but
34 /// the range of capacity function is the set of nonnegative integers.
35 /// It is not a polinomial time algorithm for counting the minimum cost
36 /// maximal flow, since it counts the minimum cost flow for every value 0..M
37 /// where \c M is the value of the maximal flow.
39 ///\author Attila Bernath
40 template <typename Graph, typename LengthMap, typename SupplyDemandMap>
43 typedef typename LengthMap::ValueType Length;
46 typedef typename SupplyDemandMap::ValueType SupplyDemand;
48 typedef typename Graph::Node Node;
49 typedef typename Graph::NodeIt NodeIt;
50 typedef typename Graph::Edge Edge;
51 typedef typename Graph::OutEdgeIt OutEdgeIt;
52 typedef typename Graph::template EdgeMap<int> EdgeIntMap;
54 // typedef ConstMap<Edge,int> ConstMap;
56 typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
57 typedef typename ResGraphType::Edge ResGraphEdge;
60 //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
61 typedef typename Graph::template NodeMap<Length> NodeMap;
62 const ResGraphType& res_graph;
63 // const EdgeIntMap& rev;
67 typedef typename LengthMap::KeyType KeyType;
68 typedef typename LengthMap::ValueType ValueType;
70 ValueType operator[](typename ResGraphType::Edge e) const {
71 if (res_graph.forward(e))
72 return ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);
74 return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);
77 ModLengthMap(const ResGraphType& _res_graph,
78 const LengthMap &o, const NodeMap &p) :
79 res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
87 const LengthMap& length;
88 const SupplyDemandMap& supply_demand;//supply or demand of nodes
95 //To store the potentila (dual variables)
96 typename Graph::template NodeMap<Length> potential;
97 //To store excess-deficit values
98 SupplyDemandMap excess_deficit;
107 MinCostFlow(Graph& _graph, LengthMap& _length, SupplyDemandMap& _supply_demand) : graph(_graph),
108 length(_length), supply_demand(_supply_demand), flow(_graph), potential(_graph){ }
111 ///Runs the algorithm.
113 ///Runs the algorithm.
115 ///\todo May be it does make sense to be able to start with a nonzero
116 /// feasible primal-dual solution pair as well.
119 //Resetting variables from previous runs
122 typedef typename Graph::template NodeMap<int> HeapMap;
123 typedef Heap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
124 std::greater<SupplyDemand> > HeapType;
126 //A heap for the excess nodes
127 HeapMap excess_nodes_map(graph,-1);
128 HeapType excess_nodes(excess_nodes_map);
130 //A heap for the deficit nodes
131 HeapMap deficit_nodes_map(graph,-1);
132 HeapType deficit_nodes(deficit_nodes_map);
134 //A container to store nonabundant arcs
135 list<Edge> nonabundant_arcs;
138 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
140 nonabundant_arcs.push_back(e);
143 //Initial value for delta
144 SupplyDemand delta = 0;
146 typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
148 //A union-find structure to store the abundant components
149 UFE::MapType abund_comp_map(graph);
150 UFE abundant_components(abund_comp_map);
154 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
155 excess_deficit.set(n,supply_demand[n]);
157 if (excess_deficit[n] > 0){
158 excess_nodes.push(n,excess_deficit[n]);
161 if (excess_deficit[n] < 0){
162 deficit_nodes.push(n, - excess_deficit[n]);
164 //Finding out starting value of delta
165 if (delta < abs(excess_deficit[n])){
166 delta = abs(excess_deficit[n]);
168 //Initialize the copy of the Dijkstra potential to zero
170 //Every single point is an abundant component initially
171 abundant_components.insert(n);
174 //It'll be allright as an initial value, though this value
175 //can be the maximum deficit here
176 SupplyDemand max_excess = delta;
178 //ConstMap<Edge,SupplyDemand> ConstEdgeMap;
180 //We need a residual graph which is uncapacitated
181 ResGraphType res_graph(graph, flow);
183 //An EdgeMap to tell which arcs are abundant
184 template typename Graph::EdgeMap<bool> abundant_arcs(graph);
186 //Let's construct the sugraph consisting only of the abundant edges
187 typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
188 ConstNodeMap const_true_map(true);
189 typedef SubGraphWrapper< Graph, ConstNodeMap,
190 template typename Graph::EdgeMap<bool> >
192 AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
194 //Let's construct the residual graph for the abundant graph
195 typedef ResGraphWrapper<const AbundantGraph,int,CapacityMap,EdgeIntMap>
197 //Again uncapacitated
198 ResAbGraph res_ab_graph(abundant_graph, flow);
200 //We need things for the bfs
201 typename ResAbGraph::NodeMap<bool> bfs_reached(res_ab_graph);
202 typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>
203 bfs_pred(res_ab_graph);
204 NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy(res_ab_graph);
205 //We want to run bfs-es (more) on this graph 'res_ab_graph'
207 typename ResAbGraph::NodeMap<bool>,
208 typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>,
209 NullMap<typename ResAbGraph::Node, int> >
210 bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
212 ModLengthMap mod_length(res_graph, length, potential);
214 Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
217 while (max_excess > 0){
219 //Reset delta if still too big
220 if (8*number_of_nodes*max_excess <= delta){
226 * Beginning of the delta scaling phase
230 SupplyDemand buf=8*number_of_nodes*delta;
231 list<Edge>::iterator i = nonabundant_arcs.begin();
232 while ( i != nonabundant_arcs.end() ){
234 Node a = abundant_components.find(res_graph.head(i));
235 Node b = abundant_components.find(res_graph.tail(i));
238 abundant_components.join(a,b);
239 //We want to push the smaller
240 //Which has greater absolut value excess/deficit
241 Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
243 Node non_root = ( a == root ) ? b : a ;
244 abundant_components.makeRep(root);
245 SupplyDemand qty_to_augment = abs(excess_deficit[non_root]);
246 //Push the positive value
247 if (excess_deficit[non_root] < 0)
248 swap(root, non_root);
249 //If the non_root node has excess/deficit at all
250 if (qty_to_augment>0){
251 //Find path and augment
253 //root should be reached
255 //Augmenting on the found path
260 n = res_graph.tail(e);
261 res_graph.augment(e,qty_to_augment);
264 //We know that non_root had positive excess
265 excess_nodes[non_root] -= qty_to_augment;
266 //But what about root node
267 //It might have been positive and so became larger
268 if (excess_deficit[root]>0){
269 excess_nodes[root] += qty_to_augment;
272 //Or negative but not turned into positive
273 deficit_nodes[root] -= qty_to_augment;
276 //Update the excess_deficit map
277 excess_deficit[non_root] -= qty_to_augment;
278 excess_deficit[root] += qty_to_augment;
284 //Marci and Zsolt says I shouldn't do such things
285 nonabundant_arcs.erase(i++);
286 abundant_arcs[i] = true;
294 Node s = excess_nodes.top();
295 SupplyDemand max_excess = excess_nodes[s];
296 Node t = deficit_nodes.top();
297 if (max_excess < deficit_nodes[t]){
298 max_excess = deficit_nodes[t];
302 while(max_excess > (n-1)*delta/n){
310 /*We know from theory that t can be reached
311 if (!dijkstra.reached(t)){
312 //There are no k paths from s to t
317 //We have to change the potential
318 FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
319 potential[n] += dijkstra.distMap()[n];
323 //Augmenting on the sortest path
327 e = dijkstra.pred(n);
328 n = dijkstra.predNode(n);
329 res_graph.augment(e,delta);
331 //Let's update the total length
332 if (res_graph.forward(e))
333 total_length += length[e];
335 total_length -= length[e];
339 //Update the excess_deficit map
340 excess_deficit[s] -= delta;
341 excess_deficit[t] += delta;
344 //Update the excess_nodes heap
345 if (delta >= excess_nodes[s]){
346 if (delta > excess_nodes[s])
347 deficit_nodes.push(s,delta - excess_nodes[s]);
352 excess_nodes[s] -= delta;
354 //Update the deficit_nodes heap
355 if (delta >= deficit_nodes[t]){
356 if (delta > deficit_nodes[t])
357 excess_nodes.push(t,delta - deficit_nodes[t]);
362 deficit_nodes[t] -= delta;
364 //Dijkstra part ends here
366 //Choose s and t again
367 s = excess_nodes.top();
368 max_excess = excess_nodes[s];
369 t = deficit_nodes.top();
370 if (max_excess < deficit_nodes[t]){
371 max_excess = deficit_nodes[t];
377 * End of the delta scaling phase
380 //Whatever this means
383 /*This is not necessary here
384 //Update the max_excess
386 FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
387 if (max_excess < excess_deficit[n]){
388 max_excess = excess_deficit[n];
394 }//while(max_excess > 0)
403 ///This function gives back the total length of the found paths.
404 ///Assumes that \c run() has been run and nothing changed since then.
405 Length totalLength(){
409 ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
410 ///be called before using this function.
411 const EdgeIntMap &getFlow() const { return flow;}
413 ///Returns a const reference to the NodeMap \c potential (the dual solution).
414 /// \pre \ref run() must be called before using this function.
415 const EdgeIntMap &getPotential() const { return potential;}
417 ///This function checks, whether the given solution is optimal
418 ///Running after a \c run() should return with true
419 ///In this "state of the art" this only check optimality, doesn't bother with feasibility
421 ///\todo Is this OK here?
422 bool checkComplementarySlackness(){
425 FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
427 mod_pot = length[e]-potential[graph.head(e)]+potential[graph.tail(e)];
429 // std::cout << fl_e << std::endl;
430 if (0<fl_e && fl_e<capacity[e]){
435 if (mod_pot > 0 && fl_e != 0)
437 if (mod_pot < 0 && fl_e != capacity[e])
445 }; //class MinCostFlow
451 #endif //HUGO_MINCOSTFLOW_H