95 .capacityMap(upper) |
101 .capacityMap(upper) |
96 .boundMaps(lower, upper) |
102 .boundMaps(lower, upper) |
97 .costMap(cost) |
103 .costMap(cost) |
98 .supplyMap(sup) |
104 .supplyMap(sup) |
99 .stSupply(n, n, k) |
105 .stSupply(n, n, k) |
|
106 .flowMap(flow) |
|
107 .potentialMap(pot) |
100 .run(); |
108 .run(); |
101 |
109 |
102 const typename MCF::FlowMap &fm = mcf.flowMap(); |
110 const MCF& const_mcf = mcf; |
103 const typename MCF::PotentialMap &pm = mcf.potentialMap(); |
111 |
104 |
112 const typename MCF::FlowMap &fm = const_mcf.flowMap(); |
105 v = mcf.totalCost(); |
113 const typename MCF::PotentialMap &pm = const_mcf.potentialMap(); |
106 double x = mcf.template totalCost<double>(); |
114 |
107 v = mcf.flow(a); |
115 v = const_mcf.totalCost(); |
108 v = mcf.potential(n); |
116 double x = const_mcf.template totalCost<double>(); |
109 mcf.flowMap(flow); |
117 v = const_mcf.flow(a); |
110 mcf.potentialMap(pot); |
118 v = const_mcf.potential(n); |
111 |
119 |
112 ignore_unused_variable_warning(fm); |
120 ignore_unused_variable_warning(fm); |
113 ignore_unused_variable_warning(pm); |
121 ignore_unused_variable_warning(pm); |
114 ignore_unused_variable_warning(x); |
122 ignore_unused_variable_warning(x); |
115 } |
123 } |
140 |
148 |
141 // Check the feasibility of the given flow (primal soluiton) |
149 // Check the feasibility of the given flow (primal soluiton) |
142 template < typename GR, typename LM, typename UM, |
150 template < typename GR, typename LM, typename UM, |
143 typename SM, typename FM > |
151 typename SM, typename FM > |
144 bool checkFlow( const GR& gr, const LM& lower, const UM& upper, |
152 bool checkFlow( const GR& gr, const LM& lower, const UM& upper, |
145 const SM& supply, const FM& flow ) |
153 const SM& supply, const FM& flow, |
|
154 ProblemType type = EQ ) |
146 { |
155 { |
147 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
156 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
148 |
157 |
149 for (ArcIt e(gr); e != INVALID; ++e) { |
158 for (ArcIt e(gr); e != INVALID; ++e) { |
150 if (flow[e] < lower[e] || flow[e] > upper[e]) return false; |
159 if (flow[e] < lower[e] || flow[e] > upper[e]) return false; |
154 typename SM::Value sum = 0; |
163 typename SM::Value sum = 0; |
155 for (OutArcIt e(gr, n); e != INVALID; ++e) |
164 for (OutArcIt e(gr, n); e != INVALID; ++e) |
156 sum += flow[e]; |
165 sum += flow[e]; |
157 for (InArcIt e(gr, n); e != INVALID; ++e) |
166 for (InArcIt e(gr, n); e != INVALID; ++e) |
158 sum -= flow[e]; |
167 sum -= flow[e]; |
159 if (sum != supply[n]) return false; |
168 bool b = (type == EQ && sum == supply[n]) || |
|
169 (type == GEQ && sum >= supply[n]) || |
|
170 (type == LEQ && sum <= supply[n]); |
|
171 if (!b) return false; |
160 } |
172 } |
161 |
173 |
162 return true; |
174 return true; |
163 } |
175 } |
164 |
176 |
165 // Check the feasibility of the given potentials (dual soluiton) |
177 // Check the feasibility of the given potentials (dual soluiton) |
166 // using the "Complementary Slackness" optimality condition |
178 // using the "Complementary Slackness" optimality condition |
167 template < typename GR, typename LM, typename UM, |
179 template < typename GR, typename LM, typename UM, |
168 typename CM, typename FM, typename PM > |
180 typename CM, typename SM, typename FM, typename PM > |
169 bool checkPotential( const GR& gr, const LM& lower, const UM& upper, |
181 bool checkPotential( const GR& gr, const LM& lower, const UM& upper, |
170 const CM& cost, const FM& flow, const PM& pi ) |
182 const CM& cost, const SM& supply, const FM& flow, |
|
183 const PM& pi ) |
171 { |
184 { |
172 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
185 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
173 |
186 |
174 bool opt = true; |
187 bool opt = true; |
175 for (ArcIt e(gr); opt && e != INVALID; ++e) { |
188 for (ArcIt e(gr); opt && e != INVALID; ++e) { |
177 cost[e] + pi[gr.source(e)] - pi[gr.target(e)]; |
190 cost[e] + pi[gr.source(e)] - pi[gr.target(e)]; |
178 opt = red_cost == 0 || |
191 opt = red_cost == 0 || |
179 (red_cost > 0 && flow[e] == lower[e]) || |
192 (red_cost > 0 && flow[e] == lower[e]) || |
180 (red_cost < 0 && flow[e] == upper[e]); |
193 (red_cost < 0 && flow[e] == upper[e]); |
181 } |
194 } |
|
195 |
|
196 for (NodeIt n(gr); opt && n != INVALID; ++n) { |
|
197 typename SM::Value sum = 0; |
|
198 for (OutArcIt e(gr, n); e != INVALID; ++e) |
|
199 sum += flow[e]; |
|
200 for (InArcIt e(gr, n); e != INVALID; ++e) |
|
201 sum -= flow[e]; |
|
202 opt = (sum == supply[n]) || (pi[n] == 0); |
|
203 } |
|
204 |
182 return opt; |
205 return opt; |
183 } |
206 } |
184 |
207 |
185 // Run a minimum cost flow algorithm and check the results |
208 // Run a minimum cost flow algorithm and check the results |
186 template < typename MCF, typename GR, |
209 template < typename MCF, typename GR, |
188 typename CM, typename SM > |
211 typename CM, typename SM > |
189 void checkMcf( const MCF& mcf, bool mcf_result, |
212 void checkMcf( const MCF& mcf, bool mcf_result, |
190 const GR& gr, const LM& lower, const UM& upper, |
213 const GR& gr, const LM& lower, const UM& upper, |
191 const CM& cost, const SM& supply, |
214 const CM& cost, const SM& supply, |
192 bool result, typename CM::Value total, |
215 bool result, typename CM::Value total, |
193 const std::string &test_id = "" ) |
216 const std::string &test_id = "", |
|
217 ProblemType type = EQ ) |
194 { |
218 { |
195 check(mcf_result == result, "Wrong result " + test_id); |
219 check(mcf_result == result, "Wrong result " + test_id); |
196 if (result) { |
220 if (result) { |
197 check(checkFlow(gr, lower, upper, supply, mcf.flowMap()), |
221 check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type), |
198 "The flow is not feasible " + test_id); |
222 "The flow is not feasible " + test_id); |
199 check(mcf.totalCost() == total, "The flow is not optimal " + test_id); |
223 check(mcf.totalCost() == total, "The flow is not optimal " + test_id); |
200 check(checkPotential(gr, lower, upper, cost, mcf.flowMap(), |
224 check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(), |
201 mcf.potentialMap()), |
225 mcf.potentialMap()), |
202 "Wrong potentials " + test_id); |
226 "Wrong potentials " + test_id); |
203 } |
227 } |
204 } |
228 } |
205 |
229 |
224 DIGRAPH_TYPEDEFS(ListDigraph); |
248 DIGRAPH_TYPEDEFS(ListDigraph); |
225 |
249 |
226 // Read the test digraph |
250 // Read the test digraph |
227 Digraph gr; |
251 Digraph gr; |
228 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr); |
252 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr); |
229 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr); |
253 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr); |
230 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); |
254 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); |
231 Node v, w; |
255 Node v, w; |
232 |
256 |
233 std::istringstream input(test_lgf); |
257 std::istringstream input(test_lgf); |
234 DigraphReader<Digraph>(gr, input) |
258 DigraphReader<Digraph>(gr, input) |
237 .arcMap("low1", l1) |
261 .arcMap("low1", l1) |
238 .arcMap("low2", l2) |
262 .arcMap("low2", l2) |
239 .nodeMap("sup1", s1) |
263 .nodeMap("sup1", s1) |
240 .nodeMap("sup2", s2) |
264 .nodeMap("sup2", s2) |
241 .nodeMap("sup3", s3) |
265 .nodeMap("sup3", s3) |
|
266 .nodeMap("sup4", s4) |
|
267 .nodeMap("sup5", s5) |
242 .node("source", v) |
268 .node("source", v) |
243 .node("target", w) |
269 .node("target", w) |
244 .run(); |
270 .run(); |
245 |
271 |
246 // A. Test NetworkSimplex with the default pivot rule |
272 // A. Test NetworkSimplex with the default pivot rule |
247 { |
273 { |
248 NetworkSimplex<Digraph> mcf(gr); |
274 NetworkSimplex<Digraph> mcf(gr); |
249 |
275 |
|
276 // Check the equality form |
250 mcf.upperMap(u).costMap(c); |
277 mcf.upperMap(u).costMap(c); |
251 checkMcf(mcf, mcf.supplyMap(s1).run(), |
278 checkMcf(mcf, mcf.supplyMap(s1).run(), |
252 gr, l1, u, c, s1, true, 5240, "#A1"); |
279 gr, l1, u, c, s1, true, 5240, "#A1"); |
253 checkMcf(mcf, mcf.stSupply(v, w, 27).run(), |
280 checkMcf(mcf, mcf.stSupply(v, w, 27).run(), |
254 gr, l1, u, c, s2, true, 7620, "#A2"); |
281 gr, l1, u, c, s2, true, 7620, "#A2"); |
265 mcf.reset(); |
292 mcf.reset(); |
266 checkMcf(mcf, mcf.run(), |
293 checkMcf(mcf, mcf.run(), |
267 gr, l1, cu, cc, s3, true, 0, "#A7"); |
294 gr, l1, cu, cc, s3, true, 0, "#A7"); |
268 checkMcf(mcf, mcf.boundMaps(l2, u).run(), |
295 checkMcf(mcf, mcf.boundMaps(l2, u).run(), |
269 gr, l2, u, cc, s3, false, 0, "#A8"); |
296 gr, l2, u, cc, s3, false, 0, "#A8"); |
|
297 |
|
298 // Check the GEQ form |
|
299 mcf.reset().upperMap(u).costMap(c).supplyMap(s4); |
|
300 checkMcf(mcf, mcf.run(), |
|
301 gr, l1, u, c, s4, true, 3530, "#A9", GEQ); |
|
302 mcf.problemType(mcf.GEQ); |
|
303 checkMcf(mcf, mcf.lowerMap(l2).run(), |
|
304 gr, l2, u, c, s4, true, 4540, "#A10", GEQ); |
|
305 mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5); |
|
306 checkMcf(mcf, mcf.run(), |
|
307 gr, l2, u, c, s5, false, 0, "#A11", GEQ); |
|
308 |
|
309 // Check the LEQ form |
|
310 mcf.reset().problemType(mcf.LEQ); |
|
311 mcf.upperMap(u).costMap(c).supplyMap(s5); |
|
312 checkMcf(mcf, mcf.run(), |
|
313 gr, l1, u, c, s5, true, 5080, "#A12", LEQ); |
|
314 checkMcf(mcf, mcf.lowerMap(l2).run(), |
|
315 gr, l2, u, c, s5, true, 5930, "#A13", LEQ); |
|
316 mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4); |
|
317 checkMcf(mcf, mcf.run(), |
|
318 gr, l2, u, c, s4, false, 0, "#A14", LEQ); |
270 } |
319 } |
271 |
320 |
272 // B. Test NetworkSimplex with each pivot rule |
321 // B. Test NetworkSimplex with each pivot rule |
273 { |
322 { |
274 NetworkSimplex<Digraph> mcf(gr); |
323 NetworkSimplex<Digraph> mcf(gr); |