Mercuriallemon-main / comparison
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | comparison | raw | help
test/min_cost_flow_test.cc
changeset 665 e652b6f9a29f
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);
lemon-main