1 /* -*- C++ -*- |
|
2 * src/lemon/min_cost_flow.h - Part of LEMON, a generic C++ optimization library |
|
3 * |
|
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
5 * (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
6 * |
|
7 * Permission to use, modify and distribute this software is granted |
|
8 * provided that this copyright notice appears in all copies. For |
|
9 * precise terms see the accompanying LICENSE file. |
|
10 * |
|
11 * This software is provided "AS IS" with no warranty of any kind, |
|
12 * express or implied, and with no claim as to its suitability for any |
|
13 * purpose. |
|
14 * |
|
15 */ |
|
16 |
|
17 #ifndef LEMON_MIN_COST_FLOW_H |
|
18 #define LEMON_MIN_COST_FLOW_H |
|
19 |
|
20 ///\ingroup flowalgs |
|
21 ///\file |
|
22 ///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost |
|
23 |
|
24 |
|
25 #include <lemon/dijkstra.h> |
|
26 #include <lemon/graph_adaptor.h> |
|
27 #include <lemon/maps.h> |
|
28 #include <vector> |
|
29 |
|
30 namespace lemon { |
|
31 |
|
32 /// \addtogroup flowalgs |
|
33 /// @{ |
|
34 |
|
35 ///\brief Implementation of an algorithm for finding a flow of value \c k |
|
36 ///(for small values of \c k) having minimal total cost between 2 nodes |
|
37 /// |
|
38 /// |
|
39 /// The class \ref lemon::MinCostFlow "MinCostFlow" implements an |
|
40 /// algorithm for finding a flow of value \c k having minimal total |
|
41 /// cost from a given source node to a given target node in an |
|
42 /// edge-weighted directed graph. To this end, the edge-capacities |
|
43 /// and edge-weights have to be nonnegative. The edge-capacities |
|
44 /// should be integers, but the edge-weights can be integers, reals |
|
45 /// or of other comparable numeric type. This algorithm is intended |
|
46 /// to be used only for small values of \c k, since it is only |
|
47 /// polynomial in k, not in the length of k (which is log k): in |
|
48 /// order to find the minimum cost flow of value \c k it finds the |
|
49 /// minimum cost flow of value \c i for every \c i between 0 and \c |
|
50 /// k. |
|
51 /// |
|
52 ///\param Graph The directed graph type the algorithm runs on. |
|
53 ///\param LengthMap The type of the length map. |
|
54 ///\param CapacityMap The capacity map type. |
|
55 /// |
|
56 ///\author Attila Bernath |
|
57 template <typename Graph, typename LengthMap, typename CapacityMap> |
|
58 class MinCostFlow { |
|
59 |
|
60 typedef typename LengthMap::Value Length; |
|
61 |
|
62 //Warning: this should be integer type |
|
63 typedef typename CapacityMap::Value Capacity; |
|
64 |
|
65 typedef typename Graph::Node Node; |
|
66 typedef typename Graph::NodeIt NodeIt; |
|
67 typedef typename Graph::Edge Edge; |
|
68 typedef typename Graph::OutEdgeIt OutEdgeIt; |
|
69 typedef typename Graph::template EdgeMap<int> EdgeIntMap; |
|
70 |
|
71 typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW; |
|
72 typedef typename ResGW::Edge ResGraphEdge; |
|
73 |
|
74 protected: |
|
75 |
|
76 const Graph& g; |
|
77 const LengthMap& length; |
|
78 const CapacityMap& capacity; |
|
79 |
|
80 EdgeIntMap flow; |
|
81 typedef typename Graph::template NodeMap<Length> PotentialMap; |
|
82 PotentialMap potential; |
|
83 |
|
84 Node s; |
|
85 Node t; |
|
86 |
|
87 Length total_length; |
|
88 |
|
89 class ModLengthMap { |
|
90 typedef typename Graph::template NodeMap<Length> NodeMap; |
|
91 const ResGW& g; |
|
92 const LengthMap &length; |
|
93 const NodeMap &pot; |
|
94 public : |
|
95 typedef typename LengthMap::Key Key; |
|
96 typedef typename LengthMap::Value Value; |
|
97 |
|
98 ModLengthMap(const ResGW& _g, |
|
99 const LengthMap &_length, const NodeMap &_pot) : |
|
100 g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { } |
|
101 |
|
102 Value operator[](typename ResGW::Edge e) const { |
|
103 if (g.forward(e)) |
|
104 return length[e]-(pot[g.target(e)]-pot[g.source(e)]); |
|
105 else |
|
106 return -length[e]-(pot[g.target(e)]-pot[g.source(e)]); |
|
107 } |
|
108 |
|
109 }; //ModLengthMap |
|
110 |
|
111 ResGW res_graph; |
|
112 ModLengthMap mod_length; |
|
113 Dijkstra<ResGW, ModLengthMap> dijkstra; |
|
114 |
|
115 public : |
|
116 |
|
117 /*! \brief The constructor of the class. |
|
118 |
|
119 \param _g The directed graph the algorithm runs on. |
|
120 \param _length The length (weight or cost) of the edges. |
|
121 \param _cap The capacity of the edges. |
|
122 \param _s Source node. |
|
123 \param _t Target node. |
|
124 */ |
|
125 MinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap, |
|
126 Node _s, Node _t) : |
|
127 g(_g), length(_length), capacity(_cap), flow(_g), potential(_g), |
|
128 s(_s), t(_t), |
|
129 res_graph(g, capacity, flow), |
|
130 mod_length(res_graph, length, potential), |
|
131 dijkstra(res_graph, mod_length) { |
|
132 reset(); |
|
133 } |
|
134 |
|
135 /*! Tries to augment the flow between s and t by 1. |
|
136 The return value shows if the augmentation is successful. |
|
137 */ |
|
138 bool augment() { |
|
139 dijkstra.run(s); |
|
140 if (!dijkstra.reached(t)) { |
|
141 |
|
142 //Unsuccessful augmentation. |
|
143 return false; |
|
144 } else { |
|
145 |
|
146 //We have to change the potential |
|
147 for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n) |
|
148 potential.set(n, potential[n]+dijkstra.distMap()[n]); |
|
149 |
|
150 //Augmenting on the shortest path |
|
151 Node n=t; |
|
152 ResGraphEdge e; |
|
153 while (n!=s){ |
|
154 e = dijkstra.pred(n); |
|
155 n = dijkstra.predNode(n); |
|
156 res_graph.augment(e,1); |
|
157 //Let's update the total length |
|
158 if (res_graph.forward(e)) |
|
159 total_length += length[e]; |
|
160 else |
|
161 total_length -= length[e]; |
|
162 } |
|
163 |
|
164 return true; |
|
165 } |
|
166 } |
|
167 |
|
168 /*! \brief Runs the algorithm. |
|
169 |
|
170 Runs the algorithm. |
|
171 Returns k if there is a flow of value at least k from s to t. |
|
172 Otherwise it returns the maximum value of a flow from s to t. |
|
173 |
|
174 \param k The value of the flow we are looking for. |
|
175 |
|
176 \todo May be it does make sense to be able to start with a nonzero |
|
177 feasible primal-dual solution pair as well. |
|
178 |
|
179 \todo If the actual flow value is bigger than k, then everything is |
|
180 cleared and the algorithm starts from zero flow. Is it a good approach? |
|
181 */ |
|
182 int run(int k) { |
|
183 if (flowValue()>k) reset(); |
|
184 while (flowValue()<k && augment()) { } |
|
185 return flowValue(); |
|
186 } |
|
187 |
|
188 /*! \brief The class is reset to zero flow and potential. |
|
189 The class is reset to zero flow and potential. |
|
190 */ |
|
191 void reset() { |
|
192 total_length=0; |
|
193 for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0); |
|
194 for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0); |
|
195 } |
|
196 |
|
197 /*! Returns the value of the actual flow. |
|
198 */ |
|
199 int flowValue() const { |
|
200 int i=0; |
|
201 for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e]; |
|
202 for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e]; |
|
203 return i; |
|
204 } |
|
205 |
|
206 /// Total weight of the found flow. |
|
207 |
|
208 /// This function gives back the total weight of the found flow. |
|
209 Length totalLength(){ |
|
210 return total_length; |
|
211 } |
|
212 |
|
213 ///Returns a const reference to the EdgeMap \c flow. |
|
214 |
|
215 ///Returns a const reference to the EdgeMap \c flow. |
|
216 const EdgeIntMap &getFlow() const { return flow;} |
|
217 |
|
218 /*! \brief Returns a const reference to the NodeMap \c potential (the dual solution). |
|
219 |
|
220 Returns a const reference to the NodeMap \c potential (the dual solution). |
|
221 */ |
|
222 const PotentialMap &getPotential() const { return potential;} |
|
223 |
|
224 /*! \brief Checking the complementary slackness optimality criteria. |
|
225 |
|
226 This function checks, whether the given flow and potential |
|
227 satisfy the complementary slackness conditions (i.e. these are optimal). |
|
228 This function only checks optimality, doesn't bother with feasibility. |
|
229 For testing purpose. |
|
230 */ |
|
231 bool checkComplementarySlackness(){ |
|
232 Length mod_pot; |
|
233 Length fl_e; |
|
234 for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
|
235 //C^{\Pi}_{i,j} |
|
236 mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)]; |
|
237 fl_e = flow[e]; |
|
238 if (0<fl_e && fl_e<capacity[e]) { |
|
239 /// \todo better comparison is needed for real types, moreover, |
|
240 /// this comparison here is superfluous. |
|
241 if (mod_pot != 0) |
|
242 return false; |
|
243 } |
|
244 else { |
|
245 if (mod_pot > 0 && fl_e != 0) |
|
246 return false; |
|
247 if (mod_pot < 0 && fl_e != capacity[e]) |
|
248 return false; |
|
249 } |
|
250 } |
|
251 return true; |
|
252 } |
|
253 |
|
254 }; //class MinCostFlow |
|
255 |
|
256 ///@} |
|
257 |
|
258 } //namespace lemon |
|
259 |
|
260 #endif //LEMON_MIN_COST_FLOW_H |
|