test/min_cost_flow_test.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 13 Nov 2009 00:23:07 +0100
changeset 884 bc75ee2ad082
parent 716 4faca85d40e6
child 885 d93490b861e9
permissions -rw-r--r--
Rework the MCF test file to help extending it (#180)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #include <iostream>
    20 #include <fstream>
    21 #include <limits>
    22 
    23 #include <lemon/list_graph.h>
    24 #include <lemon/lgf_reader.h>
    25 
    26 #include <lemon/network_simplex.h>
    27 
    28 #include <lemon/concepts/digraph.h>
    29 #include <lemon/concept_check.h>
    30 
    31 #include "test_tools.h"
    32 
    33 using namespace lemon;
    34 
    35 // Test networks
    36 char test_lgf[] =
    37   "@nodes\n"
    38   "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
    39   "    1    20   27    0   30   20   30\n"
    40   "    2    -4    0    0    0   -8   -3\n"
    41   "    3     0    0    0    0    0    0\n"
    42   "    4     0    0    0    0    0    0\n"
    43   "    5     9    0    0    0    6   11\n"
    44   "    6    -6    0    0    0   -5   -6\n"
    45   "    7     0    0    0    0    0    0\n"
    46   "    8     0    0    0    0    0    3\n"
    47   "    9     3    0    0    0    0    0\n"
    48   "   10    -2    0    0    0   -7   -2\n"
    49   "   11     0    0    0    0  -10    0\n"
    50   "   12   -20  -27    0  -30  -30  -20\n"
    51   "\n"                
    52   "@arcs\n"
    53   "       cost  cap low1 low2 low3\n"
    54   " 1  2    70   11    0    8    8\n"
    55   " 1  3   150    3    0    1    0\n"
    56   " 1  4    80   15    0    2    2\n"
    57   " 2  8    80   12    0    0    0\n"
    58   " 3  5   140    5    0    3    1\n"
    59   " 4  6    60   10    0    1    0\n"
    60   " 4  7    80    2    0    0    0\n"
    61   " 4  8   110    3    0    0    0\n"
    62   " 5  7    60   14    0    0    0\n"
    63   " 5 11   120   12    0    0    0\n"
    64   " 6  3     0    3    0    0    0\n"
    65   " 6  9   140    4    0    0    0\n"
    66   " 6 10    90    8    0    0    0\n"
    67   " 7  1    30    5    0    0   -5\n"
    68   " 8 12    60   16    0    4    3\n"
    69   " 9 12    50    6    0    0    0\n"
    70   "10 12    70   13    0    5    2\n"
    71   "10  2   100    7    0    0    0\n"
    72   "10  7    60   10    0    0   -3\n"
    73   "11 10    20   14    0    6  -20\n"
    74   "12 11    30   10    0    0  -10\n"
    75   "\n"
    76   "@attributes\n"
    77   "source 1\n"
    78   "target 12\n";
    79 
    80 char test_neg1_lgf[] =
    81   "@nodes\n"
    82   "label   sup\n"
    83   "    1   100\n"
    84   "    2     0\n"
    85   "    3     0\n"
    86   "    4  -100\n"
    87   "    5     0\n"
    88   "    6     0\n"
    89   "    7     0\n"
    90   "@arcs\n"
    91   "      cost   low1   low2\n"
    92   "1 2    100      0      0\n"
    93   "1 3     30      0      0\n"
    94   "2 4     20      0      0\n"
    95   "3 4     80      0      0\n"
    96   "3 2     50      0      0\n"
    97   "5 3     10      0      0\n"
    98   "5 6     80      0   1000\n"
    99   "6 7     30      0  -1000\n"
   100   "7 5   -120      0      0\n";
   101   
   102 char test_neg2_lgf[] =
   103   "@nodes\n"
   104   "label   sup\n"
   105   "    1   100\n"
   106   "    2  -300\n"
   107   "@arcs\n"
   108   "      cost\n"
   109   "1 2     -1\n";
   110 
   111 
   112 // Test data
   113 typedef ListDigraph Digraph;
   114 DIGRAPH_TYPEDEFS(ListDigraph);
   115 
   116 Digraph gr;
   117 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
   118 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
   119 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   120 Node v, w;
   121 
   122 Digraph neg1_gr;
   123 Digraph::ArcMap<int> neg1_c(neg1_gr), neg1_l1(neg1_gr), neg1_l2(neg1_gr);
   124 ConstMap<Arc, int> neg1_u1(std::numeric_limits<int>::max()), neg1_u2(5000);
   125 Digraph::NodeMap<int> neg1_s(neg1_gr);
   126 
   127 Digraph neg2_gr;
   128 Digraph::ArcMap<int> neg2_c(neg2_gr);
   129 ConstMap<Arc, int> neg2_l(0), neg2_u(1000);
   130 Digraph::NodeMap<int> neg2_s(neg2_gr);
   131 
   132 
   133 enum SupplyType {
   134   EQ,
   135   GEQ,
   136   LEQ
   137 };
   138 
   139 
   140 // Check the interface of an MCF algorithm
   141 template <typename GR, typename Value, typename Cost>
   142 class McfClassConcept
   143 {
   144 public:
   145 
   146   template <typename MCF>
   147   struct Constraints {
   148     void constraints() {
   149       checkConcept<concepts::Digraph, GR>();
   150       
   151       const Constraints& me = *this;
   152 
   153       MCF mcf(me.g);
   154       const MCF& const_mcf = mcf;
   155 
   156       b = mcf.reset()
   157              .lowerMap(me.lower)
   158              .upperMap(me.upper)
   159              .costMap(me.cost)
   160              .supplyMap(me.sup)
   161              .stSupply(me.n, me.n, me.k)
   162              .run();
   163 
   164       c = const_mcf.totalCost();
   165       x = const_mcf.template totalCost<double>();
   166       v = const_mcf.flow(me.a);
   167       c = const_mcf.potential(me.n);
   168       const_mcf.flowMap(fm);
   169       const_mcf.potentialMap(pm);
   170     }
   171 
   172     typedef typename GR::Node Node;
   173     typedef typename GR::Arc Arc;
   174     typedef concepts::ReadMap<Node, Value> NM;
   175     typedef concepts::ReadMap<Arc, Value> VAM;
   176     typedef concepts::ReadMap<Arc, Cost> CAM;
   177     typedef concepts::WriteMap<Arc, Value> FlowMap;
   178     typedef concepts::WriteMap<Node, Cost> PotMap;
   179   
   180     GR g;
   181     VAM lower;
   182     VAM upper;
   183     CAM cost;
   184     NM sup;
   185     Node n;
   186     Arc a;
   187     Value k;
   188 
   189     FlowMap fm;
   190     PotMap pm;
   191     bool b;
   192     double x;
   193     typename MCF::Value v;
   194     typename MCF::Cost c;
   195   };
   196 
   197 };
   198 
   199 
   200 // Check the feasibility of the given flow (primal soluiton)
   201 template < typename GR, typename LM, typename UM,
   202            typename SM, typename FM >
   203 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
   204                 const SM& supply, const FM& flow,
   205                 SupplyType type = EQ )
   206 {
   207   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   208 
   209   for (ArcIt e(gr); e != INVALID; ++e) {
   210     if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
   211   }
   212 
   213   for (NodeIt n(gr); n != INVALID; ++n) {
   214     typename SM::Value sum = 0;
   215     for (OutArcIt e(gr, n); e != INVALID; ++e)
   216       sum += flow[e];
   217     for (InArcIt e(gr, n); e != INVALID; ++e)
   218       sum -= flow[e];
   219     bool b = (type ==  EQ && sum == supply[n]) ||
   220              (type == GEQ && sum >= supply[n]) ||
   221              (type == LEQ && sum <= supply[n]);
   222     if (!b) return false;
   223   }
   224 
   225   return true;
   226 }
   227 
   228 // Check the feasibility of the given potentials (dual soluiton)
   229 // using the "Complementary Slackness" optimality condition
   230 template < typename GR, typename LM, typename UM,
   231            typename CM, typename SM, typename FM, typename PM >
   232 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   233                      const CM& cost, const SM& supply, const FM& flow, 
   234                      const PM& pi, SupplyType type )
   235 {
   236   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   237 
   238   bool opt = true;
   239   for (ArcIt e(gr); opt && e != INVALID; ++e) {
   240     typename CM::Value red_cost =
   241       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
   242     opt = red_cost == 0 ||
   243           (red_cost > 0 && flow[e] == lower[e]) ||
   244           (red_cost < 0 && flow[e] == upper[e]);
   245   }
   246   
   247   for (NodeIt n(gr); opt && n != INVALID; ++n) {
   248     typename SM::Value sum = 0;
   249     for (OutArcIt e(gr, n); e != INVALID; ++e)
   250       sum += flow[e];
   251     for (InArcIt e(gr, n); e != INVALID; ++e)
   252       sum -= flow[e];
   253     if (type != LEQ) {
   254       opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
   255     } else {
   256       opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
   257     }
   258   }
   259   
   260   return opt;
   261 }
   262 
   263 // Check whether the dual cost is equal to the primal cost
   264 template < typename GR, typename LM, typename UM,
   265            typename CM, typename SM, typename PM >
   266 bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
   267                     const CM& cost, const SM& supply, const PM& pi,
   268                     typename CM::Value total )
   269 {
   270   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   271 
   272   typename CM::Value dual_cost = 0;
   273   SM red_supply(gr);
   274   for (NodeIt n(gr); n != INVALID; ++n) {
   275     red_supply[n] = supply[n];
   276   }
   277   for (ArcIt a(gr); a != INVALID; ++a) {
   278     if (lower[a] != 0) {
   279       dual_cost += lower[a] * cost[a];
   280       red_supply[gr.source(a)] -= lower[a];
   281       red_supply[gr.target(a)] += lower[a];
   282     }
   283   }
   284   
   285   for (NodeIt n(gr); n != INVALID; ++n) {
   286     dual_cost -= red_supply[n] * pi[n];
   287   }
   288   for (ArcIt a(gr); a != INVALID; ++a) {
   289     typename CM::Value red_cost =
   290       cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
   291     dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
   292   }
   293   
   294   return dual_cost == total;
   295 }
   296 
   297 // Run a minimum cost flow algorithm and check the results
   298 template < typename MCF, typename GR,
   299            typename LM, typename UM,
   300            typename CM, typename SM,
   301            typename PT >
   302 void checkMcf( const MCF& mcf, PT mcf_result,
   303                const GR& gr, const LM& lower, const UM& upper,
   304                const CM& cost, const SM& supply,
   305                PT result, bool optimal, typename CM::Value total,
   306                const std::string &test_id = "",
   307                SupplyType type = EQ )
   308 {
   309   check(mcf_result == result, "Wrong result " + test_id);
   310   if (optimal) {
   311     typename GR::template ArcMap<typename SM::Value> flow(gr);
   312     typename GR::template NodeMap<typename CM::Value> pi(gr);
   313     mcf.flowMap(flow);
   314     mcf.potentialMap(pi);
   315     check(checkFlow(gr, lower, upper, supply, flow, type),
   316           "The flow is not feasible " + test_id);
   317     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   318     check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
   319           "Wrong potentials " + test_id);
   320     check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
   321           "Wrong dual cost " + test_id);
   322   }
   323 }
   324 
   325 template < typename MCF, typename Param >
   326 void runMcfGeqTests( Param param,
   327                      const std::string &test_str = "",
   328                      bool full_neg_cost_support = false )
   329 {
   330   MCF mcf1(gr), mcf2(neg1_gr), mcf3(neg2_gr);
   331   
   332   // Basic tests
   333   mcf1.upperMap(u).costMap(c).supplyMap(s1);
   334   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s1,
   335            mcf1.OPTIMAL, true,     5240, test_str + "-1");
   336   mcf1.stSupply(v, w, 27);
   337   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s2,
   338            mcf1.OPTIMAL, true,     7620, test_str + "-2");
   339   mcf1.lowerMap(l2).supplyMap(s1);
   340   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s1,
   341            mcf1.OPTIMAL, true,     5970, test_str + "-3");
   342   mcf1.stSupply(v, w, 27);
   343   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s2,
   344            mcf1.OPTIMAL, true,     8010, test_str + "-4");
   345   mcf1.reset().supplyMap(s1);
   346   checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s1,
   347            mcf1.OPTIMAL, true,       74, test_str + "-5");
   348   mcf1.lowerMap(l2).stSupply(v, w, 27);
   349   checkMcf(mcf1, mcf1.run(param), gr, l2, cu, cc, s2,
   350            mcf1.OPTIMAL, true,       94, test_str + "-6");
   351   mcf1.reset();
   352   checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s3,
   353            mcf1.OPTIMAL, true,        0, test_str + "-7");
   354   mcf1.lowerMap(l2).upperMap(u);
   355   checkMcf(mcf1, mcf1.run(param), gr, l2, u, cc, s3,
   356            mcf1.INFEASIBLE, false,    0, test_str + "-8");
   357   mcf1.lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
   358   checkMcf(mcf1, mcf1.run(param), gr, l3, u, c, s4,
   359            mcf1.OPTIMAL, true,     6360, test_str + "-9");
   360 
   361   // Tests for the GEQ form
   362   mcf1.reset().upperMap(u).costMap(c).supplyMap(s5);
   363   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s5,
   364            mcf1.OPTIMAL, true,     3530, test_str + "-10", GEQ);
   365   mcf1.lowerMap(l2);
   366   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
   367            mcf1.OPTIMAL, true,     4540, test_str + "-11", GEQ);
   368   mcf1.supplyMap(s6);
   369   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
   370            mcf1.INFEASIBLE, false,    0, test_str + "-12", GEQ);
   371 
   372   // Tests with negative costs
   373   mcf2.lowerMap(neg1_l1).costMap(neg1_c).supplyMap(neg1_s);
   374   checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u1, neg1_c, neg1_s,
   375            mcf2.UNBOUNDED, false,     0, test_str + "-13");
   376   mcf2.upperMap(neg1_u2);
   377   checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u2, neg1_c, neg1_s,
   378            mcf2.OPTIMAL, true,   -40000, test_str + "-14");
   379   mcf2.reset().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s);
   380   checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l2, neg1_u1, neg1_c, neg1_s,
   381            mcf2.UNBOUNDED, false,     0, test_str + "-15");
   382 
   383   mcf3.costMap(neg2_c).supplyMap(neg2_s);
   384   if (full_neg_cost_support) {
   385     checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
   386              mcf3.OPTIMAL, true,   -300, test_str + "-16", GEQ);
   387   } else {
   388     checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
   389              mcf3.UNBOUNDED, false,   0, test_str + "-17", GEQ);
   390   }
   391   mcf3.upperMap(neg2_u);
   392   checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
   393            mcf3.OPTIMAL, true,     -300, test_str + "-18", GEQ);
   394 }
   395 
   396 template < typename MCF, typename Param >
   397 void runMcfLeqTests( Param param,
   398                      const std::string &test_str = "" )
   399 {
   400   // Tests for the LEQ form
   401   MCF mcf1(gr);
   402   mcf1.supplyType(mcf1.LEQ);
   403   mcf1.upperMap(u).costMap(c).supplyMap(s6);
   404   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6,
   405            mcf1.OPTIMAL, true,   5080, test_str + "-19", LEQ);
   406   mcf1.lowerMap(l2);
   407   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
   408            mcf1.OPTIMAL, true,   5930, test_str + "-20", LEQ);
   409   mcf1.supplyMap(s5);
   410   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
   411            mcf1.INFEASIBLE, false,  0, test_str + "-21", LEQ);
   412 }
   413 
   414 
   415 int main()
   416 {
   417   // Read the test networks
   418   std::istringstream input(test_lgf);
   419   DigraphReader<Digraph>(gr, input)
   420     .arcMap("cost", c)
   421     .arcMap("cap", u)
   422     .arcMap("low1", l1)
   423     .arcMap("low2", l2)
   424     .arcMap("low3", l3)
   425     .nodeMap("sup1", s1)
   426     .nodeMap("sup2", s2)
   427     .nodeMap("sup3", s3)
   428     .nodeMap("sup4", s4)
   429     .nodeMap("sup5", s5)
   430     .nodeMap("sup6", s6)
   431     .node("source", v)
   432     .node("target", w)
   433     .run();
   434   
   435   std::istringstream neg_inp1(test_neg1_lgf);
   436   DigraphReader<Digraph>(neg1_gr, neg_inp1)
   437     .arcMap("cost", neg1_c)
   438     .arcMap("low1", neg1_l1)
   439     .arcMap("low2", neg1_l2)
   440     .nodeMap("sup", neg1_s)
   441     .run();
   442   
   443   std::istringstream neg_inp2(test_neg2_lgf);
   444   DigraphReader<Digraph>(neg2_gr, neg_inp2)
   445     .arcMap("cost", neg2_c)
   446     .nodeMap("sup", neg2_s)
   447     .run();
   448   
   449   // Check the interface of NetworkSimplex
   450   {
   451     typedef concepts::Digraph GR;
   452     checkConcept< McfClassConcept<GR, int, int>,
   453                   NetworkSimplex<GR> >();
   454     checkConcept< McfClassConcept<GR, double, double>,
   455                   NetworkSimplex<GR, double> >();
   456     checkConcept< McfClassConcept<GR, int, double>,
   457                   NetworkSimplex<GR, int, double> >();
   458   }
   459 
   460   // Test NetworkSimplex
   461   { 
   462     typedef NetworkSimplex<Digraph> MCF;
   463     runMcfGeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE", true);
   464     runMcfLeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE");
   465     runMcfGeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE", true);
   466     runMcfLeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE");
   467     runMcfGeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS", true);
   468     runMcfLeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS");
   469     runMcfGeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL", true);
   470     runMcfLeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL");
   471     runMcfGeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL", true);
   472     runMcfLeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL");
   473   }
   474 
   475   return 0;
   476 }