# source:lemon-0.x/src/work/athos/mincostflow.h@635:933f593824c2

Last change on this file since 635:933f593824c2 was 635:933f593824c2, checked in by athos, 20 years ago

Started mincostflow.

File size: 8.1 KB
Line
1// -*- c++ -*-
2#ifndef HUGO_MINCOSTFLOW_H
3#define HUGO_MINCOSTFLOW_H
4
5///\ingroup galgs
6///\file
7///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
8
9
10#include <hugo/dijkstra.h>
11#include <hugo/graph_wrapper.h>
12#include <hugo/maps.h>
13#include <vector>
14#include <for_each_macros.h>
15
16namespace hugo {
17
19/// @{
20
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
23  ///
24  ///
25  /// The class \ref hugo::MinCostFlow "MinCostFlow" implements
26  /// an algorithm for solving the following general minimum cost flow problem>
27  ///
28  ///
29  ///
30  /// \warning It is assumed here that the problem has a feasible solution
31  ///
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.
37  ///
38  ///\author Attila Bernath
39  template <typename Graph, typename LengthMap, typename SupplyDemandMap>
40  class MinCostFlow {
41
42    typedef typename LengthMap::ValueType Length;
43
44
45    typedef typename SupplyDemandMap::ValueType SupplyDemand;
46
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;
52
53    //    typedef ConstMap<Edge,int> ConstMap;
54
55    typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
56    typedef typename ResGraphType::Edge ResGraphEdge;
57
58    class ModLengthMap {
59      //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
60      typedef typename Graph::template NodeMap<Length> NodeMap;
61      const ResGraphType& G;
62      //      const EdgeIntMap& rev;
63      const LengthMap &ol;
64      const NodeMap &pot;
65    public :
66      typedef typename LengthMap::KeyType KeyType;
67      typedef typename LengthMap::ValueType ValueType;
68
69      ValueType operator[](typename ResGraphType::Edge e) const {
70        if (G.forward(e))
72        else
74      }
75
76      ModLengthMap(const ResGraphType& _G,
77                   const LengthMap &o,  const NodeMap &p) :
78        G(_G), /*rev(_rev),*/ ol(o), pot(p){};
79    };//ModLengthMap
80
81
82  protected:
83
84    //Input
85    const Graph& G;
86    const LengthMap& length;
87    const SupplyDemandMap& supply_demand;//supply or demand of nodes
88
89
90    //auxiliary variables
91
92    //To store the flow
93    EdgeIntMap flow;
94    //To store the potentila (dual variables)
95    typename Graph::template NodeMap<Length> potential;
96    //To store excess-deficit values
97    SupplyDemandMap excess_deficit;
98
99
100    Length total_length;
101
102
103  public :
104
105
106    MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G),
107      length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
108
109
110    ///Runs the algorithm.
111
112    ///Runs the algorithm.
113
114    ///\todo May be it does make sense to be able to start with a nonzero
115    /// feasible primal-dual solution pair as well.
116    int run() {
117
118      //Resetting variables from previous runs
119      //total_length = 0;
120
121      typedef typename Graph::template NodeMap<int> HeapMap;
122      typedef Heap<Node, SupplyDemand, typename Graph::template NodeMap<int>,
123        std::greater<SupplyDemand> >    HeapType;
124
125      //A heap for the excess nodes
126      HeapMap excess_nodes_map(G,-1);
127      HeapType excess_nodes(excess_nodes_map);
128
129      //A heap for the deficit nodes
130      HeapMap deficit_nodes_map(G,-1);
131      HeapType deficit_nodes(deficit_nodes_map);
132
133
134      FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
135        flow.set(e,0);
136      }
137
138      //Initial value for delta
139      SupplyDemand delta = 0;
140
141      FOR_EACH_LOC(typename Graph::NodeIt, n, G){
142        excess_deficit.set(n,supply_demand[n]);
143        //A supply node
144        if (excess_deficit[n] > 0){
145          excess_nodes.push(n,excess_deficit[n]);
146        }
147        //A demand node
148        if (excess_deficit[n] < 0){
149          deficit_nodes.push(n, - excess_deficit[n]);
150        }
151        //Finding out starting value of delta
152        if (delta < abs(excess_deficit[n])){
153          delta = abs(excess_deficit[n]);
154        }
155        //Initialize the copy of the Dijkstra potential to zero
156        potential.set(n,0);
157      }
158
159      //It'll be allright as an initial value, though this value
160      //can be the maximum deficit here
161      SupplyDemand max_excess = delta;
162
163      //We need a residual graph which is uncapacitated
164      ResGraphType res_graph(G, flow);
165
166
167
168      ModLengthMap mod_length(res_graph, length, potential);
169
170      Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
171
172
173      while (max_excess > 0){
174
175
176        //Merge and stuff
177
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];
183        }
184
185
186        while(max_excess > ){
187
188
189          //s es t valasztasa
190
191          //Dijkstra part
192          dijkstra.run(s);
193
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
197            break;
198          };
199          */
200
201          //We have to change the potential
202          FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
203            potential[n] += dijkstra.distMap()[n];
204          }
205
206
207          //Augmenting on the sortest path
208          Node n=t;
209          ResGraphEdge e;
210          while (n!=s){
211            e = dijkstra.pred(n);
212            n = dijkstra.predNode(n);
213            res_graph.augment(e,delta);
214            /*
215            //Let's update the total length
216            if (res_graph.forward(e))
217              total_length += length[e];
218            else
219              total_length -= length[e];
220            */
221          }
222
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]);
227            excess_nodes.pop();
228
229          }
230          else{
231            excess_nodes[s] -= delta;
232          }
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]);
237            deficit_nodes.pop();
238
239          }
240          else{
241            deficit_nodes[t] -= delta;
242          }
243          //Dijkstra part ends here
244        }
245
246        /*
247         * End of the delta scaling phase
248        */
249
250        //Whatever this means
251        delta = delta / 2;
252
253        /*This is not necessary here
254        //Update the max_excess
255        max_excess = 0;
256        FOR_EACH_LOC(typename Graph::NodeIt, n, G){
257          if (max_excess < excess_deficit[n]){
258            max_excess = excess_deficit[n];
259          }
260        }
261        */
262        //Reset delta if still too big
263        if (8*number_of_nodes*max_excess <= delta){
264          delta = max_excess;
265
266        }
267
268      }//while(max_excess > 0)
269
270
271      return i;
272    }
273
274
275
276
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(){
281    }
282
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;}
286
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;}
290
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
294    ///
295    ///\todo Is this OK here?
296    bool checkComplementarySlackness(){
297      Length mod_pot;
298      Length fl_e;
299      FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
300        //C^{\Pi}_{i,j}
302        fl_e = flow[e];
303        //      std::cout << fl_e << std::endl;
304        if (0<fl_e && fl_e<capacity[e]){
305          if (mod_pot != 0)
306            return false;
307        }
308        else{
309          if (mod_pot > 0 && fl_e != 0)
310            return false;
311          if (mod_pot < 0 && fl_e != capacity[e])
312            return false;
313        }
314      }
315      return true;
316    }
317
318
319  }; //class MinCostFlow
320
321  ///@}
322
323} //namespace hugo
324
325#endif //HUGO_MINCOSTFLOW_H
Note: See TracBrowser for help on using the repository browser.