Started mincostflow.
2 #ifndef HUGO_MINCOSTFLOW_H
3 #define HUGO_MINCOSTFLOW_H
7 ///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
10 #include <hugo/dijkstra.h>
11 #include <hugo/graph_wrapper.h>
12 #include <hugo/maps.h>
14 #include <for_each_macros.h>
21 ///\brief Implementation of an algorithm for finding a flow of value \c k
22 ///(for small values of \c k) having minimal total cost between 2 nodes
25 /// The class \ref hugo::MinCostFlow "MinCostFlow" implements
26 /// an algorithm for solving the following general minimum cost flow problem>
30 /// \warning It is assumed here that the problem has a feasible solution
32 /// The range of the length (weight) function is nonnegative reals but
33 /// the range of capacity function is the set of nonnegative integers.
34 /// It is not a polinomial time algorithm for counting the minimum cost
35 /// maximal flow, since it counts the minimum cost flow for every value 0..M
36 /// where \c M is the value of the maximal flow.
38 ///\author Attila Bernath
39 template <typename Graph, typename LengthMap, typename SupplyDemandMap>
42 typedef typename LengthMap::ValueType Length;
45 typedef typename SupplyDemandMap::ValueType SupplyDemand;
47 typedef typename Graph::Node Node;
48 typedef typename Graph::NodeIt NodeIt;
49 typedef typename Graph::Edge Edge;
50 typedef typename Graph::OutEdgeIt OutEdgeIt;
51 typedef typename Graph::template EdgeMap<int> EdgeIntMap;
53 // typedef ConstMap<Edge,int> ConstMap;
55 typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
56 typedef typename ResGraphType::Edge ResGraphEdge;
59 //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
60 typedef typename Graph::template NodeMap<Length> NodeMap;
61 const ResGraphType& G;
62 // const EdgeIntMap& rev;
66 typedef typename LengthMap::KeyType KeyType;
67 typedef typename LengthMap::ValueType ValueType;
69 ValueType operator[](typename ResGraphType::Edge e) const {
71 return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
73 return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
76 ModLengthMap(const ResGraphType& _G,
77 const LengthMap &o, const NodeMap &p) :
78 G(_G), /*rev(_rev),*/ ol(o), pot(p){};
86 const LengthMap& length;
87 const SupplyDemandMap& supply_demand;//supply or demand of nodes
94 //To store the potentila (dual variables)
95 typename Graph::template NodeMap<Length> potential;
96 //To store excess-deficit values
97 SupplyDemandMap excess_deficit;
106 MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G),
107 length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
110 ///Runs the algorithm.
112 ///Runs the algorithm.
114 ///\todo May be it does make sense to be able to start with a nonzero
115 /// feasible primal-dual solution pair as well.
118 //Resetting variables from previous runs
121 typedef typename Graph::template NodeMap<int> HeapMap;
122 typedef Heap<Node, SupplyDemand, typename Graph::template NodeMap<int>,
123 std::greater<SupplyDemand> > HeapType;
125 //A heap for the excess nodes
126 HeapMap excess_nodes_map(G,-1);
127 HeapType excess_nodes(excess_nodes_map);
129 //A heap for the deficit nodes
130 HeapMap deficit_nodes_map(G,-1);
131 HeapType deficit_nodes(deficit_nodes_map);
134 FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
138 //Initial value for delta
139 SupplyDemand delta = 0;
141 FOR_EACH_LOC(typename Graph::NodeIt, n, G){
142 excess_deficit.set(n,supply_demand[n]);
144 if (excess_deficit[n] > 0){
145 excess_nodes.push(n,excess_deficit[n]);
148 if (excess_deficit[n] < 0){
149 deficit_nodes.push(n, - excess_deficit[n]);
151 //Finding out starting value of delta
152 if (delta < abs(excess_deficit[n])){
153 delta = abs(excess_deficit[n]);
155 //Initialize the copy of the Dijkstra potential to zero
159 //It'll be allright as an initial value, though this value
160 //can be the maximum deficit here
161 SupplyDemand max_excess = delta;
163 //We need a residual graph which is uncapacitated
164 ResGraphType res_graph(G, flow);
168 ModLengthMap mod_length(res_graph, length, potential);
170 Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
173 while (max_excess > 0){
178 Node s = excess_nodes.top();
179 SupplyDemand max_excess = excess_nodes[s];
180 Node t = deficit_nodes.top();
181 if (max_excess < dificit_nodes[t]){
182 max_excess = dificit_nodes[t];
186 while(max_excess > ){
194 /*We know from theory that t can be reached
195 if (!dijkstra.reached(t)){
196 //There are no k paths from s to t
201 //We have to change the potential
202 FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
203 potential[n] += dijkstra.distMap()[n];
207 //Augmenting on the sortest path
211 e = dijkstra.pred(n);
212 n = dijkstra.predNode(n);
213 res_graph.augment(e,delta);
215 //Let's update the total length
216 if (res_graph.forward(e))
217 total_length += length[e];
219 total_length -= length[e];
223 //Update the excess_nodes heap
224 if (delta >= excess_nodes[s]){
225 if (delta > excess_nodes[s])
226 deficit_nodes.push(s,delta - excess_nodes[s]);
231 excess_nodes[s] -= delta;
233 //Update the deficit_nodes heap
234 if (delta >= deficit_nodes[t]){
235 if (delta > deficit_nodes[t])
236 excess_nodes.push(t,delta - deficit_nodes[t]);
241 deficit_nodes[t] -= delta;
243 //Dijkstra part ends here
247 * End of the delta scaling phase
250 //Whatever this means
253 /*This is not necessary here
254 //Update the max_excess
256 FOR_EACH_LOC(typename Graph::NodeIt, n, G){
257 if (max_excess < excess_deficit[n]){
258 max_excess = excess_deficit[n];
262 //Reset delta if still too big
263 if (8*number_of_nodes*max_excess <= delta){
268 }//while(max_excess > 0)
277 ///This function gives back the total length of the found paths.
278 ///Assumes that \c run() has been run and nothing changed since then.
279 Length totalLength(){
283 ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
284 ///be called before using this function.
285 const EdgeIntMap &getFlow() const { return flow;}
287 ///Returns a const reference to the NodeMap \c potential (the dual solution).
288 /// \pre \ref run() must be called before using this function.
289 const EdgeIntMap &getPotential() const { return potential;}
291 ///This function checks, whether the given solution is optimal
292 ///Running after a \c run() should return with true
293 ///In this "state of the art" this only check optimality, doesn't bother with feasibility
295 ///\todo Is this OK here?
296 bool checkComplementarySlackness(){
299 FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
301 mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
303 // std::cout << fl_e << std::endl;
304 if (0<fl_e && fl_e<capacity[e]){
309 if (mod_pot > 0 && fl_e != 0)
311 if (mod_pot < 0 && fl_e != capacity[e])
319 }; //class MinCostFlow
325 #endif //HUGO_MINCOSTFLOW_H