r711:cc61d09f053b
Tue, 12 May 2009 12:08:06
[Not Reviewed]
0 comments (0 inline)
0
1
0
| ... | ... |
@@ -165,72 +165,112 @@ |
| 165 | 165 |
if (!b) return false; |
| 166 | 166 |
} |
| 167 | 167 |
|
| 168 | 168 |
return true; |
| 169 | 169 |
} |
| 170 | 170 |
|
| 171 | 171 |
// Check the feasibility of the given potentials (dual soluiton) |
| 172 | 172 |
// using the "Complementary Slackness" optimality condition |
| 173 | 173 |
template < typename GR, typename LM, typename UM, |
| 174 | 174 |
typename CM, typename SM, typename FM, typename PM > |
| 175 | 175 |
bool checkPotential( const GR& gr, const LM& lower, const UM& upper, |
| 176 | 176 |
const CM& cost, const SM& supply, const FM& flow, |
| 177 |
const PM& pi ) |
|
| 177 |
const PM& pi, SupplyType type ) |
|
| 178 | 178 |
{
|
| 179 | 179 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
| 180 | 180 |
|
| 181 | 181 |
bool opt = true; |
| 182 | 182 |
for (ArcIt e(gr); opt && e != INVALID; ++e) {
|
| 183 | 183 |
typename CM::Value red_cost = |
| 184 | 184 |
cost[e] + pi[gr.source(e)] - pi[gr.target(e)]; |
| 185 | 185 |
opt = red_cost == 0 || |
| 186 | 186 |
(red_cost > 0 && flow[e] == lower[e]) || |
| 187 | 187 |
(red_cost < 0 && flow[e] == upper[e]); |
| 188 | 188 |
} |
| 189 | 189 |
|
| 190 | 190 |
for (NodeIt n(gr); opt && n != INVALID; ++n) {
|
| 191 | 191 |
typename SM::Value sum = 0; |
| 192 | 192 |
for (OutArcIt e(gr, n); e != INVALID; ++e) |
| 193 | 193 |
sum += flow[e]; |
| 194 | 194 |
for (InArcIt e(gr, n); e != INVALID; ++e) |
| 195 | 195 |
sum -= flow[e]; |
| 196 |
|
|
| 196 |
if (type != LEQ) {
|
|
| 197 |
opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0); |
|
| 198 |
} else {
|
|
| 199 |
opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0); |
|
| 200 |
} |
|
| 197 | 201 |
} |
| 198 | 202 |
|
| 199 | 203 |
return opt; |
| 200 | 204 |
} |
| 201 | 205 |
|
| 206 |
// Check whether the dual cost is equal to the primal cost |
|
| 207 |
template < typename GR, typename LM, typename UM, |
|
| 208 |
typename CM, typename SM, typename PM > |
|
| 209 |
bool checkDualCost( const GR& gr, const LM& lower, const UM& upper, |
|
| 210 |
const CM& cost, const SM& supply, const PM& pi, |
|
| 211 |
typename CM::Value total ) |
|
| 212 |
{
|
|
| 213 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
|
| 214 |
|
|
| 215 |
typename CM::Value dual_cost = 0; |
|
| 216 |
SM red_supply(gr); |
|
| 217 |
for (NodeIt n(gr); n != INVALID; ++n) {
|
|
| 218 |
red_supply[n] = supply[n]; |
|
| 219 |
} |
|
| 220 |
for (ArcIt a(gr); a != INVALID; ++a) {
|
|
| 221 |
if (lower[a] != 0) {
|
|
| 222 |
dual_cost += lower[a] * cost[a]; |
|
| 223 |
red_supply[gr.source(a)] -= lower[a]; |
|
| 224 |
red_supply[gr.target(a)] += lower[a]; |
|
| 225 |
} |
|
| 226 |
} |
|
| 227 |
|
|
| 228 |
for (NodeIt n(gr); n != INVALID; ++n) {
|
|
| 229 |
dual_cost -= red_supply[n] * pi[n]; |
|
| 230 |
} |
|
| 231 |
for (ArcIt a(gr); a != INVALID; ++a) {
|
|
| 232 |
typename CM::Value red_cost = |
|
| 233 |
cost[a] + pi[gr.source(a)] - pi[gr.target(a)]; |
|
| 234 |
dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0); |
|
| 235 |
} |
|
| 236 |
|
|
| 237 |
return dual_cost == total; |
|
| 238 |
} |
|
| 239 |
|
|
| 202 | 240 |
// Run a minimum cost flow algorithm and check the results |
| 203 | 241 |
template < typename MCF, typename GR, |
| 204 | 242 |
typename LM, typename UM, |
| 205 | 243 |
typename CM, typename SM, |
| 206 | 244 |
typename PT > |
| 207 | 245 |
void checkMcf( const MCF& mcf, PT mcf_result, |
| 208 | 246 |
const GR& gr, const LM& lower, const UM& upper, |
| 209 | 247 |
const CM& cost, const SM& supply, |
| 210 | 248 |
PT result, bool optimal, typename CM::Value total, |
| 211 | 249 |
const std::string &test_id = "", |
| 212 | 250 |
SupplyType type = EQ ) |
| 213 | 251 |
{
|
| 214 | 252 |
check(mcf_result == result, "Wrong result " + test_id); |
| 215 | 253 |
if (optimal) {
|
| 216 | 254 |
typename GR::template ArcMap<typename SM::Value> flow(gr); |
| 217 | 255 |
typename GR::template NodeMap<typename CM::Value> pi(gr); |
| 218 | 256 |
mcf.flowMap(flow); |
| 219 | 257 |
mcf.potentialMap(pi); |
| 220 | 258 |
check(checkFlow(gr, lower, upper, supply, flow, type), |
| 221 | 259 |
"The flow is not feasible " + test_id); |
| 222 | 260 |
check(mcf.totalCost() == total, "The flow is not optimal " + test_id); |
| 223 |
check(checkPotential(gr, lower, upper, cost, supply, flow, pi), |
|
| 261 |
check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type), |
|
| 224 | 262 |
"Wrong potentials " + test_id); |
| 263 |
check(checkDualCost(gr, lower, upper, cost, supply, pi, total), |
|
| 264 |
"Wrong dual cost " + test_id); |
|
| 225 | 265 |
} |
| 226 | 266 |
} |
| 227 | 267 |
|
| 228 | 268 |
int main() |
| 229 | 269 |
{
|
| 230 | 270 |
// Check the interfaces |
| 231 | 271 |
{
|
| 232 | 272 |
typedef concepts::Digraph GR; |
| 233 | 273 |
checkConcept< McfClassConcept<GR, int, int>, |
| 234 | 274 |
NetworkSimplex<GR> >(); |
| 235 | 275 |
checkConcept< McfClassConcept<GR, double, double>, |
| 236 | 276 |
NetworkSimplex<GR, double> >(); |
| ... | ... |
@@ -257,63 +297,74 @@ |
| 257 | 297 |
.arcMap("low2", l2)
|
| 258 | 298 |
.arcMap("low3", l3)
|
| 259 | 299 |
.nodeMap("sup1", s1)
|
| 260 | 300 |
.nodeMap("sup2", s2)
|
| 261 | 301 |
.nodeMap("sup3", s3)
|
| 262 | 302 |
.nodeMap("sup4", s4)
|
| 263 | 303 |
.nodeMap("sup5", s5)
|
| 264 | 304 |
.nodeMap("sup6", s6)
|
| 265 | 305 |
.node("source", v)
|
| 266 | 306 |
.node("target", w)
|
| 267 | 307 |
.run(); |
| 268 | 308 |
|
| 269 |
// Build a test digraph for testing negative costs |
|
| 270 |
Digraph ngr; |
|
| 271 |
Node n1 = ngr.addNode(); |
|
| 272 |
Node n2 = ngr.addNode(); |
|
| 273 |
Node n3 = ngr.addNode(); |
|
| 274 |
Node n4 = ngr.addNode(); |
|
| 275 |
Node n5 = ngr.addNode(); |
|
| 276 |
Node n6 = ngr.addNode(); |
|
| 277 |
|
|
| 309 |
// Build test digraphs with negative costs |
|
| 310 |
Digraph neg_gr; |
|
| 311 |
Node n1 = neg_gr.addNode(); |
|
| 312 |
Node n2 = neg_gr.addNode(); |
|
| 313 |
Node n3 = neg_gr.addNode(); |
|
| 314 |
Node n4 = neg_gr.addNode(); |
|
| 315 |
Node n5 = neg_gr.addNode(); |
|
| 316 |
Node n6 = neg_gr.addNode(); |
|
| 317 |
Node n7 = neg_gr.addNode(); |
|
| 278 | 318 |
|
| 279 |
Arc a1 = ngr.addArc(n1, n2); |
|
| 280 |
Arc a2 = ngr.addArc(n1, n3); |
|
| 281 |
Arc a3 = ngr.addArc(n2, n4); |
|
| 282 |
Arc a4 = ngr.addArc(n3, n4); |
|
| 283 |
Arc a5 = ngr.addArc(n3, n2); |
|
| 284 |
Arc a6 = ngr.addArc(n5, n3); |
|
| 285 |
Arc a7 = ngr.addArc(n5, n6); |
|
| 286 |
Arc a8 = ngr.addArc(n6, n7); |
|
| 287 |
Arc |
|
| 319 |
Arc a1 = neg_gr.addArc(n1, n2); |
|
| 320 |
Arc a2 = neg_gr.addArc(n1, n3); |
|
| 321 |
Arc a3 = neg_gr.addArc(n2, n4); |
|
| 322 |
Arc a4 = neg_gr.addArc(n3, n4); |
|
| 323 |
Arc a5 = neg_gr.addArc(n3, n2); |
|
| 324 |
Arc a6 = neg_gr.addArc(n5, n3); |
|
| 325 |
Arc a7 = neg_gr.addArc(n5, n6); |
|
| 326 |
Arc a8 = neg_gr.addArc(n6, n7); |
|
| 327 |
Arc a9 = neg_gr.addArc(n7, n5); |
|
| 288 | 328 |
|
| 289 |
Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0); |
|
| 290 |
ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000); |
|
| 291 |
Digraph:: |
|
| 329 |
Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0); |
|
| 330 |
ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000); |
|
| 331 |
Digraph::NodeMap<int> neg_s(neg_gr, 0); |
|
| 292 | 332 |
|
| 293 |
nl2[a7] = 1000; |
|
| 294 |
nl2[a8] = -1000; |
|
| 333 |
neg_l2[a7] = 1000; |
|
| 334 |
neg_l2[a8] = -1000; |
|
| 295 | 335 |
|
| 296 |
ns[n1] = 100; |
|
| 297 |
ns[n4] = -100; |
|
| 336 |
neg_s[n1] = 100; |
|
| 337 |
neg_s[n4] = -100; |
|
| 298 | 338 |
|
| 299 |
nc[a1] = 100; |
|
| 300 |
nc[a2] = 30; |
|
| 301 |
nc[a3] = 20; |
|
| 302 |
nc[a4] = 80; |
|
| 303 |
nc[a5] = 50; |
|
| 304 |
nc[a6] = 10; |
|
| 305 |
nc[a7] = 80; |
|
| 306 |
nc[a8] = 30; |
|
| 307 |
|
|
| 339 |
neg_c[a1] = 100; |
|
| 340 |
neg_c[a2] = 30; |
|
| 341 |
neg_c[a3] = 20; |
|
| 342 |
neg_c[a4] = 80; |
|
| 343 |
neg_c[a5] = 50; |
|
| 344 |
neg_c[a6] = 10; |
|
| 345 |
neg_c[a7] = 80; |
|
| 346 |
neg_c[a8] = 30; |
|
| 347 |
neg_c[a9] = -120; |
|
| 348 |
|
|
| 349 |
Digraph negs_gr; |
|
| 350 |
Digraph::NodeMap<int> negs_s(negs_gr); |
|
| 351 |
Digraph::ArcMap<int> negs_c(negs_gr); |
|
| 352 |
ConstMap<Arc, int> negs_l(0), negs_u(1000); |
|
| 353 |
n1 = negs_gr.addNode(); |
|
| 354 |
n2 = negs_gr.addNode(); |
|
| 355 |
negs_s[n1] = 100; |
|
| 356 |
negs_s[n2] = -300; |
|
| 357 |
negs_c[negs_gr.addArc(n1, n2)] = -1; |
|
| 358 |
|
|
| 308 | 359 |
|
| 309 | 360 |
// A. Test NetworkSimplex with the default pivot rule |
| 310 | 361 |
{
|
| 311 | 362 |
NetworkSimplex<Digraph> mcf(gr); |
| 312 | 363 |
|
| 313 | 364 |
// Check the equality form |
| 314 | 365 |
mcf.upperMap(u).costMap(c); |
| 315 | 366 |
checkMcf(mcf, mcf.supplyMap(s1).run(), |
| 316 | 367 |
gr, l1, u, c, s1, mcf.OPTIMAL, true, 5240, "#A1"); |
| 317 | 368 |
checkMcf(mcf, mcf.stSupply(v, w, 27).run(), |
| 318 | 369 |
gr, l1, u, c, s2, mcf.OPTIMAL, true, 7620, "#A2"); |
| 319 | 370 |
mcf.lowerMap(l2); |
| ... | ... |
@@ -333,49 +384,55 @@ |
| 333 | 384 |
gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8"); |
| 334 | 385 |
mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4); |
| 335 | 386 |
checkMcf(mcf, mcf.run(), |
| 336 | 387 |
gr, l3, u, c, s4, mcf.OPTIMAL, true, 6360, "#A9"); |
| 337 | 388 |
|
| 338 | 389 |
// Check the GEQ form |
| 339 | 390 |
mcf.reset().upperMap(u).costMap(c).supplyMap(s5); |
| 340 | 391 |
checkMcf(mcf, mcf.run(), |
| 341 | 392 |
gr, l1, u, c, s5, mcf.OPTIMAL, true, 3530, "#A10", GEQ); |
| 342 | 393 |
mcf.supplyType(mcf.GEQ); |
| 343 | 394 |
checkMcf(mcf, mcf.lowerMap(l2).run(), |
| 344 | 395 |
gr, l2, u, c, s5, mcf.OPTIMAL, true, 4540, "#A11", GEQ); |
| 345 |
mcf. |
|
| 396 |
mcf.supplyMap(s6); |
|
| 346 | 397 |
checkMcf(mcf, mcf.run(), |
| 347 | 398 |
gr, l2, u, c, s6, mcf.INFEASIBLE, false, 0, "#A12", GEQ); |
| 348 | 399 |
|
| 349 | 400 |
// Check the LEQ form |
| 350 | 401 |
mcf.reset().supplyType(mcf.LEQ); |
| 351 | 402 |
mcf.upperMap(u).costMap(c).supplyMap(s6); |
| 352 | 403 |
checkMcf(mcf, mcf.run(), |
| 353 | 404 |
gr, l1, u, c, s6, mcf.OPTIMAL, true, 5080, "#A13", LEQ); |
| 354 | 405 |
checkMcf(mcf, mcf.lowerMap(l2).run(), |
| 355 | 406 |
gr, l2, u, c, s6, mcf.OPTIMAL, true, 5930, "#A14", LEQ); |
| 356 |
mcf. |
|
| 407 |
mcf.supplyMap(s5); |
|
| 357 | 408 |
checkMcf(mcf, mcf.run(), |
| 358 | 409 |
gr, l2, u, c, s5, mcf.INFEASIBLE, false, 0, "#A15", LEQ); |
| 359 | 410 |
|
| 360 | 411 |
// Check negative costs |
| 361 |
NetworkSimplex<Digraph> nmcf(ngr); |
|
| 362 |
nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns); |
|
| 363 |
checkMcf(nmcf, nmcf.run(), |
|
| 364 |
ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A16"); |
|
| 365 |
checkMcf(nmcf, nmcf.upperMap(nu2).run(), |
|
| 366 |
ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true, -40000, "#A17"); |
|
| 367 |
nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns); |
|
| 368 |
checkMcf(nmcf, nmcf.run(), |
|
| 369 |
|
|
| 412 |
NetworkSimplex<Digraph> neg_mcf(neg_gr); |
|
| 413 |
neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s); |
|
| 414 |
checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1, |
|
| 415 |
neg_c, neg_s, neg_mcf.UNBOUNDED, false, 0, "#A16"); |
|
| 416 |
neg_mcf.upperMap(neg_u2); |
|
| 417 |
checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2, |
|
| 418 |
neg_c, neg_s, neg_mcf.OPTIMAL, true, -40000, "#A17"); |
|
| 419 |
neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s); |
|
| 420 |
checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1, |
|
| 421 |
neg_c, neg_s, neg_mcf.UNBOUNDED, false, 0, "#A18"); |
|
| 422 |
|
|
| 423 |
NetworkSimplex<Digraph> negs_mcf(negs_gr); |
|
| 424 |
negs_mcf.costMap(negs_c).supplyMap(negs_s); |
|
| 425 |
checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u, |
|
| 426 |
negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ); |
|
| 370 | 427 |
} |
| 371 | 428 |
|
| 372 | 429 |
// B. Test NetworkSimplex with each pivot rule |
| 373 | 430 |
{
|
| 374 | 431 |
NetworkSimplex<Digraph> mcf(gr); |
| 375 | 432 |
mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2); |
| 376 | 433 |
|
| 377 | 434 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE), |
| 378 | 435 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B1"); |
| 379 | 436 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE), |
| 380 | 437 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B2"); |
| 381 | 438 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH), |
0 comments (0 inline)