r711:cc61d09f053b
Tue, 12 May 2009 12:08:06
[Not Reviewed]
0 comments (0 inline)
0
1
0
| ... | ... |
@@ -171,13 +171,13 @@ |
| 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 = |
| ... | ... |
@@ -190,18 +190,56 @@ |
| 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, |
| ... | ... |
@@ -217,14 +255,16 @@ |
| 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 |
| ... | ... |
@@ -263,51 +303,62 @@ |
| 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 |
| ... | ... |
@@ -339,37 +390,43 @@ |
| 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); |
0 comments (0 inline)