test/min_cost_flow_test.cc
changeset 613 b1811c363299
parent 607 9ad8d2122b50
child 615 e3d9bff447ed
equal deleted inserted replaced
3:39c7edb88ab7 4:c052bbe1cccb
    31 
    31 
    32 using namespace lemon;
    32 using namespace lemon;
    33 
    33 
    34 char test_lgf[] =
    34 char test_lgf[] =
    35   "@nodes\n"
    35   "@nodes\n"
    36   "label  sup1 sup2 sup3\n"
    36   "label  sup1 sup2 sup3 sup4 sup5\n"
    37   "    1    20   27    0\n"
    37   "    1    20   27    0   20   30\n"
    38   "    2    -4    0    0\n"
    38   "    2    -4    0    0   -8   -3\n"
    39   "    3     0    0    0\n"
    39   "    3     0    0    0    0    0\n"
    40   "    4     0    0    0\n"
    40   "    4     0    0    0    0    0\n"
    41   "    5     9    0    0\n"
    41   "    5     9    0    0    6   11\n"
    42   "    6    -6    0    0\n"
    42   "    6    -6    0    0   -5   -6\n"
    43   "    7     0    0    0\n"
    43   "    7     0    0    0    0    0\n"
    44   "    8     0    0    0\n"
    44   "    8     0    0    0    0    3\n"
    45   "    9     3    0    0\n"
    45   "    9     3    0    0    0    0\n"
    46   "   10    -2    0    0\n"
    46   "   10    -2    0    0   -7   -2\n"
    47   "   11     0    0    0\n"
    47   "   11     0    0    0  -10    0\n"
    48   "   12   -20  -27    0\n"
    48   "   12   -20  -27    0  -30  -20\n"
    49   "\n"
    49   "\n"
    50   "@arcs\n"
    50   "@arcs\n"
    51   "       cost  cap low1 low2\n"
    51   "       cost  cap low1 low2\n"
    52   " 1  2    70   11    0    8\n"
    52   " 1  2    70   11    0    8\n"
    53   " 1  3   150    3    0    1\n"
    53   " 1  3   150    3    0    1\n"
    74   "@attributes\n"
    74   "@attributes\n"
    75   "source 1\n"
    75   "source 1\n"
    76   "target 12\n";
    76   "target 12\n";
    77 
    77 
    78 
    78 
       
    79 enum ProblemType {
       
    80   EQ,
       
    81   GEQ,
       
    82   LEQ
       
    83 };
       
    84 
    79 // Check the interface of an MCF algorithm
    85 // Check the interface of an MCF algorithm
    80 template <typename GR, typename Flow, typename Cost>
    86 template <typename GR, typename Flow, typename Cost>
    81 class McfClassConcept
    87 class McfClassConcept
    82 {
    88 {
    83 public:
    89 public:
    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);