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