34 /// It is not a polinomial time algorithm for counting the minimum cost |
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 |
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. |
36 /// where \c M is the value of the maximal flow. |
37 /// |
37 /// |
38 ///\author Attila Bernath |
38 ///\author Attila Bernath |
39 template <typename Graph, typename LengthMap, typename SupplyMap> |
39 template <typename Graph, typename LengthMap, typename SupplyDemandMap> |
40 class MinCostFlow { |
40 class MinCostFlow { |
41 |
41 |
42 typedef typename LengthMap::ValueType Length; |
42 typedef typename LengthMap::ValueType Length; |
43 |
43 |
44 |
44 |
45 typedef typename SupplyMap::ValueType Supply; |
45 typedef typename SupplyDemandMap::ValueType SupplyDemand; |
46 |
46 |
47 typedef typename Graph::Node Node; |
47 typedef typename Graph::Node Node; |
48 typedef typename Graph::NodeIt NodeIt; |
48 typedef typename Graph::NodeIt NodeIt; |
49 typedef typename Graph::Edge Edge; |
49 typedef typename Graph::Edge Edge; |
50 typedef typename Graph::OutEdgeIt OutEdgeIt; |
50 typedef typename Graph::OutEdgeIt OutEdgeIt; |
82 protected: |
82 protected: |
83 |
83 |
84 //Input |
84 //Input |
85 const Graph& G; |
85 const Graph& G; |
86 const LengthMap& length; |
86 const LengthMap& length; |
87 const SupplyMap& supply;//supply or demand of nodes |
87 const SupplyDemandMap& supply_demand;//supply or demand of nodes |
88 |
88 |
89 |
89 |
90 //auxiliary variables |
90 //auxiliary variables |
91 |
91 |
92 //To store the flow |
92 //To store the flow |
93 EdgeIntMap flow; |
93 EdgeIntMap flow; |
94 //To store the potentila (dual variables) |
94 //To store the potentila (dual variables) |
95 typename Graph::template NodeMap<Length> potential; |
95 typename Graph::template NodeMap<Length> potential; |
96 //To store excess-deficit values |
96 //To store excess-deficit values |
97 SupplyMap excess; |
97 SupplyDemandMap excess_deficit; |
98 |
98 |
99 |
99 |
100 Length total_length; |
100 Length total_length; |
101 |
101 |
102 |
102 |
103 public : |
103 public : |
104 |
104 |
105 |
105 |
106 MinCostFlow(Graph& _G, LengthMap& _length, SupplyMap& _supply) : G(_G), |
106 MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G), |
107 length(_length), supply(_supply), flow(_G), potential(_G){ } |
107 length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ } |
108 |
108 |
109 |
109 |
110 ///Runs the algorithm. |
110 ///Runs the algorithm. |
111 |
111 |
112 ///Runs the algorithm. |
112 ///Runs the algorithm. |
113 ///Returns k if there are at least k edge-disjoint paths from s to t. |
113 |
114 ///Otherwise it returns the number of found edge-disjoint paths from s to t. |
|
115 ///\todo May be it does make sense to be able to start with a nonzero |
114 ///\todo May be it does make sense to be able to start with a nonzero |
116 /// feasible primal-dual solution pair as well. |
115 /// feasible primal-dual solution pair as well. |
117 int run() { |
116 int run() { |
118 |
117 |
119 //Resetting variables from previous runs |
118 //Resetting variables from previous runs |
120 total_length = 0; |
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 |
121 |
133 |
122 FOR_EACH_LOC(typename Graph::EdgeIt, e, G){ |
134 FOR_EACH_LOC(typename Graph::EdgeIt, e, G){ |
123 flow.set(e,0); |
135 flow.set(e,0); |
124 } |
136 } |
125 |
137 |
126 //Initial value for delta |
138 //Initial value for delta |
127 Supply delta = 0; |
139 SupplyDemand delta = 0; |
128 |
140 |
129 FOR_EACH_LOC(typename Graph::NodeIt, n, G){ |
141 FOR_EACH_LOC(typename Graph::NodeIt, n, G){ |
130 if (delta < abs(supply[e])){ |
142 excess_deficit.set(n,supply_demand[n]); |
131 delta = abs(supply[e]); |
143 //A supply node |
132 } |
144 if (excess_deficit[n] > 0){ |
133 excess.set(n,supply[e]); |
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 } |
134 //Initialize the copy of the Dijkstra potential to zero |
155 //Initialize the copy of the Dijkstra potential to zero |
135 potential.set(n,0); |
156 potential.set(n,0); |
136 } |
157 } |
137 |
158 |
138 |
159 //It'll be allright as an initial value, though this value |
|
160 //can be the maximum deficit here |
|
161 SupplyDemand max_excess = delta; |
139 |
162 |
140 //We need a residual graph which is uncapacitated |
163 //We need a residual graph which is uncapacitated |
141 ResGraphType res_graph(G, flow); |
164 ResGraphType res_graph(G, flow); |
142 |
165 |
143 |
166 |
145 ModLengthMap mod_length(res_graph, length, potential); |
168 ModLengthMap mod_length(res_graph, length, potential); |
146 |
169 |
147 Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length); |
170 Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length); |
148 |
171 |
149 |
172 |
150 int i; |
173 while (max_excess > 0){ |
151 for (i=0; i<k; ++i){ |
174 |
152 dijkstra.run(s); |
|
153 if (!dijkstra.reached(t)){ |
|
154 //There are no k paths from s to t |
|
155 break; |
|
156 }; |
|
157 |
175 |
158 //We have to copy the potential |
176 //Merge and stuff |
159 FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){ |
177 |
160 potential[n] += dijkstra.distMap()[n]; |
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 |
161 } |
244 } |
162 |
245 |
163 /* |
246 /* |
164 { |
247 * End of the delta scaling phase |
165 |
|
166 typename ResGraphType::NodeIt n; |
|
167 for ( res_graph.first(n) ; res_graph.valid(n) ; res_graph.next(n) ) { |
|
168 potential[n] += dijkstra.distMap()[n]; |
|
169 } |
|
170 } |
|
171 */ |
248 */ |
172 |
249 |
173 //Augmenting on the sortest path |
250 //Whatever this means |
174 Node n=t; |
251 delta = delta / 2; |
175 ResGraphEdge e; |
252 |
176 while (n!=s){ |
253 /*This is not necessary here |
177 e = dijkstra.pred(n); |
254 //Update the max_excess |
178 n = dijkstra.predNode(n); |
255 max_excess = 0; |
179 res_graph.augment(e,delta); |
256 FOR_EACH_LOC(typename Graph::NodeIt, n, G){ |
180 //Let's update the total length |
257 if (max_excess < excess_deficit[n]){ |
181 if (res_graph.forward(e)) |
258 max_excess = excess_deficit[n]; |
182 total_length += length[e]; |
259 } |
183 else |
260 } |
184 total_length -= length[e]; |
261 */ |
185 } |
262 //Reset delta if still too big |
186 |
263 if (8*number_of_nodes*max_excess <= delta){ |
|
264 delta = max_excess; |
187 |
265 |
188 } |
266 } |
|
267 |
|
268 }//while(max_excess > 0) |
189 |
269 |
190 |
270 |
191 return i; |
271 return i; |
192 } |
272 } |
193 |
273 |