test/min_cost_flow_test.cc
changeset 666 1993af615e68
parent 642 111698359429
child 669 4faca85d40e6
equal deleted inserted replaced
7:2a8458afeadb 8:7f1e238a1cf3
   172 // using the "Complementary Slackness" optimality condition
   172 // using the "Complementary Slackness" optimality condition
   173 template < typename GR, typename LM, typename UM,
   173 template < typename GR, typename LM, typename UM,
   174            typename CM, typename SM, typename FM, typename PM >
   174            typename CM, typename SM, typename FM, typename PM >
   175 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   175 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   176                      const CM& cost, const SM& supply, const FM& flow, 
   176                      const CM& cost, const SM& supply, const FM& flow, 
   177                      const PM& pi )
   177                      const PM& pi, SupplyType type )
   178 {
   178 {
   179   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   179   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   180 
   180 
   181   bool opt = true;
   181   bool opt = true;
   182   for (ArcIt e(gr); opt && e != INVALID; ++e) {
   182   for (ArcIt e(gr); opt && e != INVALID; ++e) {
   191     typename SM::Value sum = 0;
   191     typename SM::Value sum = 0;
   192     for (OutArcIt e(gr, n); e != INVALID; ++e)
   192     for (OutArcIt e(gr, n); e != INVALID; ++e)
   193       sum += flow[e];
   193       sum += flow[e];
   194     for (InArcIt e(gr, n); e != INVALID; ++e)
   194     for (InArcIt e(gr, n); e != INVALID; ++e)
   195       sum -= flow[e];
   195       sum -= flow[e];
   196     opt = (sum == supply[n]) || (pi[n] == 0);
   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   return opt;
   203   return opt;
       
   204 }
       
   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;
   200 }
   238 }
   201 
   239 
   202 // Run a minimum cost flow algorithm and check the results
   240 // Run a minimum cost flow algorithm and check the results
   203 template < typename MCF, typename GR,
   241 template < typename MCF, typename GR,
   204            typename LM, typename UM,
   242            typename LM, typename UM,
   218     mcf.flowMap(flow);
   256     mcf.flowMap(flow);
   219     mcf.potentialMap(pi);
   257     mcf.potentialMap(pi);
   220     check(checkFlow(gr, lower, upper, supply, flow, type),
   258     check(checkFlow(gr, lower, upper, supply, flow, type),
   221           "The flow is not feasible " + test_id);
   259           "The flow is not feasible " + test_id);
   222     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   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           "Wrong potentials " + test_id);
   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 int main()
   268 int main()
   229 {
   269 {
   264     .nodeMap("sup6", s6)
   304     .nodeMap("sup6", s6)
   265     .node("source", v)
   305     .node("source", v)
   266     .node("target", w)
   306     .node("target", w)
   267     .run();
   307     .run();
   268   
   308   
   269   // Build a test digraph for testing negative costs
   309   // Build test digraphs with negative costs
   270   Digraph ngr;
   310   Digraph neg_gr;
   271   Node n1 = ngr.addNode();
   311   Node n1 = neg_gr.addNode();
   272   Node n2 = ngr.addNode();
   312   Node n2 = neg_gr.addNode();
   273   Node n3 = ngr.addNode();
   313   Node n3 = neg_gr.addNode();
   274   Node n4 = ngr.addNode();
   314   Node n4 = neg_gr.addNode();
   275   Node n5 = ngr.addNode();
   315   Node n5 = neg_gr.addNode();
   276   Node n6 = ngr.addNode();
   316   Node n6 = neg_gr.addNode();
   277   Node n7 = ngr.addNode();
   317   Node n7 = neg_gr.addNode();
   278   
   318   
   279   Arc a1 = ngr.addArc(n1, n2);
   319   Arc a1 = neg_gr.addArc(n1, n2);
   280   Arc a2 = ngr.addArc(n1, n3);
   320   Arc a2 = neg_gr.addArc(n1, n3);
   281   Arc a3 = ngr.addArc(n2, n4);
   321   Arc a3 = neg_gr.addArc(n2, n4);
   282   Arc a4 = ngr.addArc(n3, n4);
   322   Arc a4 = neg_gr.addArc(n3, n4);
   283   Arc a5 = ngr.addArc(n3, n2);
   323   Arc a5 = neg_gr.addArc(n3, n2);
   284   Arc a6 = ngr.addArc(n5, n3);
   324   Arc a6 = neg_gr.addArc(n5, n3);
   285   Arc a7 = ngr.addArc(n5, n6);
   325   Arc a7 = neg_gr.addArc(n5, n6);
   286   Arc a8 = ngr.addArc(n6, n7);
   326   Arc a8 = neg_gr.addArc(n6, n7);
   287   Arc a9 = ngr.addArc(n7, n5);
   327   Arc a9 = neg_gr.addArc(n7, n5);
   288   
   328   
   289   Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
   329   Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
   290   ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
   330   ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
   291   Digraph::NodeMap<int> ns(ngr, 0);
   331   Digraph::NodeMap<int> neg_s(neg_gr, 0);
   292   
   332   
   293   nl2[a7] =  1000;
   333   neg_l2[a7] =  1000;
   294   nl2[a8] = -1000;
   334   neg_l2[a8] = -1000;
   295   
   335   
   296   ns[n1] =  100;
   336   neg_s[n1] =  100;
   297   ns[n4] = -100;
   337   neg_s[n4] = -100;
   298   
   338   
   299   nc[a1] =  100;
   339   neg_c[a1] =  100;
   300   nc[a2] =   30;
   340   neg_c[a2] =   30;
   301   nc[a3] =   20;
   341   neg_c[a3] =   20;
   302   nc[a4] =   80;
   342   neg_c[a4] =   80;
   303   nc[a5] =   50;
   343   neg_c[a5] =   50;
   304   nc[a6] =   10;
   344   neg_c[a6] =   10;
   305   nc[a7] =   80;
   345   neg_c[a7] =   80;
   306   nc[a8] =   30;
   346   neg_c[a8] =   30;
   307   nc[a9] = -120;
   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   // A. Test NetworkSimplex with the default pivot rule
   360   // A. Test NetworkSimplex with the default pivot rule
   310   {
   361   {
   311     NetworkSimplex<Digraph> mcf(gr);
   362     NetworkSimplex<Digraph> mcf(gr);
   312 
   363 
   340     checkMcf(mcf, mcf.run(),
   391     checkMcf(mcf, mcf.run(),
   341              gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
   392              gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
   342     mcf.supplyType(mcf.GEQ);
   393     mcf.supplyType(mcf.GEQ);
   343     checkMcf(mcf, mcf.lowerMap(l2).run(),
   394     checkMcf(mcf, mcf.lowerMap(l2).run(),
   344              gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   395              gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   345     mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
   396     mcf.supplyMap(s6);
   346     checkMcf(mcf, mcf.run(),
   397     checkMcf(mcf, mcf.run(),
   347              gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   398              gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   348 
   399 
   349     // Check the LEQ form
   400     // Check the LEQ form
   350     mcf.reset().supplyType(mcf.LEQ);
   401     mcf.reset().supplyType(mcf.LEQ);
   351     mcf.upperMap(u).costMap(c).supplyMap(s6);
   402     mcf.upperMap(u).costMap(c).supplyMap(s6);
   352     checkMcf(mcf, mcf.run(),
   403     checkMcf(mcf, mcf.run(),
   353              gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   404              gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   354     checkMcf(mcf, mcf.lowerMap(l2).run(),
   405     checkMcf(mcf, mcf.lowerMap(l2).run(),
   355              gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   406              gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   356     mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
   407     mcf.supplyMap(s5);
   357     checkMcf(mcf, mcf.run(),
   408     checkMcf(mcf, mcf.run(),
   358              gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   409              gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   359 
   410 
   360     // Check negative costs
   411     // Check negative costs
   361     NetworkSimplex<Digraph> nmcf(ngr);
   412     NetworkSimplex<Digraph> neg_mcf(neg_gr);
   362     nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
   413     neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
   363     checkMcf(nmcf, nmcf.run(),
   414     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
   364       ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
   415       neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
   365     checkMcf(nmcf, nmcf.upperMap(nu2).run(),
   416     neg_mcf.upperMap(neg_u2);
   366       ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
   417     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
   367     nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
   418       neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
   368     checkMcf(nmcf, nmcf.run(),
   419     neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
   369       ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
   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   // B. Test NetworkSimplex with each pivot rule
   429   // B. Test NetworkSimplex with each pivot rule
   373   {
   430   {
   374     NetworkSimplex<Digraph> mcf(gr);
   431     NetworkSimplex<Digraph> mcf(gr);