test/min_cost_flow_test.cc
author Alpar Juttner <alpar@cs.elte.hu>
Mon, 08 Aug 2011 12:36:16 +0200
branch1.1
changeset 1081 f1398882a928
parent 716 4faca85d40e6
permissions -rw-r--r--
Unify sources
     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-2011
     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 char test_lgf[] =
    36   "@nodes\n"
    37   "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
    38   "    1    20   27    0   30   20   30\n"
    39   "    2    -4    0    0    0   -8   -3\n"
    40   "    3     0    0    0    0    0    0\n"
    41   "    4     0    0    0    0    0    0\n"
    42   "    5     9    0    0    0    6   11\n"
    43   "    6    -6    0    0    0   -5   -6\n"
    44   "    7     0    0    0    0    0    0\n"
    45   "    8     0    0    0    0    0    3\n"
    46   "    9     3    0    0    0    0    0\n"
    47   "   10    -2    0    0    0   -7   -2\n"
    48   "   11     0    0    0    0  -10    0\n"
    49   "   12   -20  -27    0  -30  -30  -20\n"
    50   "\n"
    51   "@arcs\n"
    52   "       cost  cap low1 low2 low3\n"
    53   " 1  2    70   11    0    8    8\n"
    54   " 1  3   150    3    0    1    0\n"
    55   " 1  4    80   15    0    2    2\n"
    56   " 2  8    80   12    0    0    0\n"
    57   " 3  5   140    5    0    3    1\n"
    58   " 4  6    60   10    0    1    0\n"
    59   " 4  7    80    2    0    0    0\n"
    60   " 4  8   110    3    0    0    0\n"
    61   " 5  7    60   14    0    0    0\n"
    62   " 5 11   120   12    0    0    0\n"
    63   " 6  3     0    3    0    0    0\n"
    64   " 6  9   140    4    0    0    0\n"
    65   " 6 10    90    8    0    0    0\n"
    66   " 7  1    30    5    0    0   -5\n"
    67   " 8 12    60   16    0    4    3\n"
    68   " 9 12    50    6    0    0    0\n"
    69   "10 12    70   13    0    5    2\n"
    70   "10  2   100    7    0    0    0\n"
    71   "10  7    60   10    0    0   -3\n"
    72   "11 10    20   14    0    6  -20\n"
    73   "12 11    30   10    0    0  -10\n"
    74   "\n"
    75   "@attributes\n"
    76   "source 1\n"
    77   "target 12\n";
    78 
    79 
    80 enum SupplyType {
    81   EQ,
    82   GEQ,
    83   LEQ
    84 };
    85 
    86 // Check the interface of an MCF algorithm
    87 template <typename GR, typename Value, typename Cost>
    88 class McfClassConcept
    89 {
    90 public:
    91 
    92   template <typename MCF>
    93   struct Constraints {
    94     void constraints() {
    95       checkConcept<concepts::Digraph, GR>();
    96 
    97       const Constraints& me = *this;
    98 
    99       MCF mcf(me.g);
   100       const MCF& const_mcf = mcf;
   101 
   102       b = mcf.reset()
   103              .lowerMap(me.lower)
   104              .upperMap(me.upper)
   105              .costMap(me.cost)
   106              .supplyMap(me.sup)
   107              .stSupply(me.n, me.n, me.k)
   108              .run();
   109 
   110       c = const_mcf.totalCost();
   111       x = const_mcf.template totalCost<double>();
   112       v = const_mcf.flow(me.a);
   113       c = const_mcf.potential(me.n);
   114       const_mcf.flowMap(fm);
   115       const_mcf.potentialMap(pm);
   116     }
   117 
   118     typedef typename GR::Node Node;
   119     typedef typename GR::Arc Arc;
   120     typedef concepts::ReadMap<Node, Value> NM;
   121     typedef concepts::ReadMap<Arc, Value> VAM;
   122     typedef concepts::ReadMap<Arc, Cost> CAM;
   123     typedef concepts::WriteMap<Arc, Value> FlowMap;
   124     typedef concepts::WriteMap<Node, Cost> PotMap;
   125 
   126     GR g;
   127     VAM lower;
   128     VAM upper;
   129     CAM cost;
   130     NM sup;
   131     Node n;
   132     Arc a;
   133     Value k;
   134 
   135     FlowMap fm;
   136     PotMap pm;
   137     bool b;
   138     double x;
   139     typename MCF::Value v;
   140     typename MCF::Cost c;
   141   };
   142 
   143 };
   144 
   145 
   146 // Check the feasibility of the given flow (primal soluiton)
   147 template < typename GR, typename LM, typename UM,
   148            typename SM, typename FM >
   149 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
   150                 const SM& supply, const FM& flow,
   151                 SupplyType type = EQ )
   152 {
   153   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   154 
   155   for (ArcIt e(gr); e != INVALID; ++e) {
   156     if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
   157   }
   158 
   159   for (NodeIt n(gr); n != INVALID; ++n) {
   160     typename SM::Value sum = 0;
   161     for (OutArcIt e(gr, n); e != INVALID; ++e)
   162       sum += flow[e];
   163     for (InArcIt e(gr, n); e != INVALID; ++e)
   164       sum -= flow[e];
   165     bool b = (type ==  EQ && sum == supply[n]) ||
   166              (type == GEQ && sum >= supply[n]) ||
   167              (type == LEQ && sum <= supply[n]);
   168     if (!b) return false;
   169   }
   170 
   171   return true;
   172 }
   173 
   174 // Check the feasibility of the given potentials (dual soluiton)
   175 // using the "Complementary Slackness" optimality condition
   176 template < typename GR, typename LM, typename UM,
   177            typename CM, typename SM, typename FM, typename PM >
   178 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   179                      const CM& cost, const SM& supply, const FM& flow,
   180                      const PM& pi, SupplyType type )
   181 {
   182   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   183 
   184   bool opt = true;
   185   for (ArcIt e(gr); opt && e != INVALID; ++e) {
   186     typename CM::Value red_cost =
   187       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
   188     opt = red_cost == 0 ||
   189           (red_cost > 0 && flow[e] == lower[e]) ||
   190           (red_cost < 0 && flow[e] == upper[e]);
   191   }
   192 
   193   for (NodeIt n(gr); opt && n != INVALID; ++n) {
   194     typename SM::Value sum = 0;
   195     for (OutArcIt e(gr, n); e != INVALID; ++e)
   196       sum += flow[e];
   197     for (InArcIt e(gr, n); e != INVALID; ++e)
   198       sum -= flow[e];
   199     if (type != LEQ) {
   200       opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
   201     } else {
   202       opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
   203     }
   204   }
   205 
   206   return opt;
   207 }
   208 
   209 // Check whether the dual cost is equal to the primal cost
   210 template < typename GR, typename LM, typename UM,
   211            typename CM, typename SM, typename PM >
   212 bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
   213                     const CM& cost, const SM& supply, const PM& pi,
   214                     typename CM::Value total )
   215 {
   216   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   217 
   218   typename CM::Value dual_cost = 0;
   219   SM red_supply(gr);
   220   for (NodeIt n(gr); n != INVALID; ++n) {
   221     red_supply[n] = supply[n];
   222   }
   223   for (ArcIt a(gr); a != INVALID; ++a) {
   224     if (lower[a] != 0) {
   225       dual_cost += lower[a] * cost[a];
   226       red_supply[gr.source(a)] -= lower[a];
   227       red_supply[gr.target(a)] += lower[a];
   228     }
   229   }
   230 
   231   for (NodeIt n(gr); n != INVALID; ++n) {
   232     dual_cost -= red_supply[n] * pi[n];
   233   }
   234   for (ArcIt a(gr); a != INVALID; ++a) {
   235     typename CM::Value red_cost =
   236       cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
   237     dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
   238   }
   239 
   240   return dual_cost == total;
   241 }
   242 
   243 // Run a minimum cost flow algorithm and check the results
   244 template < typename MCF, typename GR,
   245            typename LM, typename UM,
   246            typename CM, typename SM,
   247            typename PT >
   248 void checkMcf( const MCF& mcf, PT mcf_result,
   249                const GR& gr, const LM& lower, const UM& upper,
   250                const CM& cost, const SM& supply,
   251                PT result, bool optimal, typename CM::Value total,
   252                const std::string &test_id = "",
   253                SupplyType type = EQ )
   254 {
   255   check(mcf_result == result, "Wrong result " + test_id);
   256   if (optimal) {
   257     typename GR::template ArcMap<typename SM::Value> flow(gr);
   258     typename GR::template NodeMap<typename CM::Value> pi(gr);
   259     mcf.flowMap(flow);
   260     mcf.potentialMap(pi);
   261     check(checkFlow(gr, lower, upper, supply, flow, type),
   262           "The flow is not feasible " + test_id);
   263     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   264     check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
   265           "Wrong potentials " + test_id);
   266     check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
   267           "Wrong dual cost " + test_id);
   268   }
   269 }
   270 
   271 int main()
   272 {
   273   // Check the interfaces
   274   {
   275     typedef concepts::Digraph GR;
   276     checkConcept< McfClassConcept<GR, int, int>,
   277                   NetworkSimplex<GR> >();
   278     checkConcept< McfClassConcept<GR, double, double>,
   279                   NetworkSimplex<GR, double> >();
   280     checkConcept< McfClassConcept<GR, int, double>,
   281                   NetworkSimplex<GR, int, double> >();
   282   }
   283 
   284   // Run various MCF tests
   285   typedef ListDigraph Digraph;
   286   DIGRAPH_TYPEDEFS(ListDigraph);
   287 
   288   // Read the test digraph
   289   Digraph gr;
   290   Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
   291   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
   292   ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   293   Node v, w;
   294 
   295   std::istringstream input(test_lgf);
   296   DigraphReader<Digraph>(gr, input)
   297     .arcMap("cost", c)
   298     .arcMap("cap", u)
   299     .arcMap("low1", l1)
   300     .arcMap("low2", l2)
   301     .arcMap("low3", l3)
   302     .nodeMap("sup1", s1)
   303     .nodeMap("sup2", s2)
   304     .nodeMap("sup3", s3)
   305     .nodeMap("sup4", s4)
   306     .nodeMap("sup5", s5)
   307     .nodeMap("sup6", s6)
   308     .node("source", v)
   309     .node("target", w)
   310     .run();
   311 
   312   // Build test digraphs with negative costs
   313   Digraph neg_gr;
   314   Node n1 = neg_gr.addNode();
   315   Node n2 = neg_gr.addNode();
   316   Node n3 = neg_gr.addNode();
   317   Node n4 = neg_gr.addNode();
   318   Node n5 = neg_gr.addNode();
   319   Node n6 = neg_gr.addNode();
   320   Node n7 = neg_gr.addNode();
   321 
   322   Arc a1 = neg_gr.addArc(n1, n2);
   323   Arc a2 = neg_gr.addArc(n1, n3);
   324   Arc a3 = neg_gr.addArc(n2, n4);
   325   Arc a4 = neg_gr.addArc(n3, n4);
   326   Arc a5 = neg_gr.addArc(n3, n2);
   327   Arc a6 = neg_gr.addArc(n5, n3);
   328   Arc a7 = neg_gr.addArc(n5, n6);
   329   Arc a8 = neg_gr.addArc(n6, n7);
   330   Arc a9 = neg_gr.addArc(n7, n5);
   331 
   332   Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
   333   ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
   334   Digraph::NodeMap<int> neg_s(neg_gr, 0);
   335 
   336   neg_l2[a7] =  1000;
   337   neg_l2[a8] = -1000;
   338 
   339   neg_s[n1] =  100;
   340   neg_s[n4] = -100;
   341 
   342   neg_c[a1] =  100;
   343   neg_c[a2] =   30;
   344   neg_c[a3] =   20;
   345   neg_c[a4] =   80;
   346   neg_c[a5] =   50;
   347   neg_c[a6] =   10;
   348   neg_c[a7] =   80;
   349   neg_c[a8] =   30;
   350   neg_c[a9] = -120;
   351 
   352   Digraph negs_gr;
   353   Digraph::NodeMap<int> negs_s(negs_gr);
   354   Digraph::ArcMap<int> negs_c(negs_gr);
   355   ConstMap<Arc, int> negs_l(0), negs_u(1000);
   356   n1 = negs_gr.addNode();
   357   n2 = negs_gr.addNode();
   358   negs_s[n1] = 100;
   359   negs_s[n2] = -300;
   360   negs_c[negs_gr.addArc(n1, n2)] = -1;
   361 
   362 
   363   // A. Test NetworkSimplex with the default pivot rule
   364   {
   365     NetworkSimplex<Digraph> mcf(gr);
   366 
   367     // Check the equality form
   368     mcf.upperMap(u).costMap(c);
   369     checkMcf(mcf, mcf.supplyMap(s1).run(),
   370              gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
   371     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   372              gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
   373     mcf.lowerMap(l2);
   374     checkMcf(mcf, mcf.supplyMap(s1).run(),
   375              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
   376     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   377              gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
   378     mcf.reset();
   379     checkMcf(mcf, mcf.supplyMap(s1).run(),
   380              gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
   381     checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
   382              gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
   383     mcf.reset();
   384     checkMcf(mcf, mcf.run(),
   385              gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
   386     checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
   387              gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
   388     mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
   389     checkMcf(mcf, mcf.run(),
   390              gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
   391 
   392     // Check the GEQ form
   393     mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
   394     checkMcf(mcf, mcf.run(),
   395              gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
   396     mcf.supplyType(mcf.GEQ);
   397     checkMcf(mcf, mcf.lowerMap(l2).run(),
   398              gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   399     mcf.supplyMap(s6);
   400     checkMcf(mcf, mcf.run(),
   401              gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   402 
   403     // Check the LEQ form
   404     mcf.reset().supplyType(mcf.LEQ);
   405     mcf.upperMap(u).costMap(c).supplyMap(s6);
   406     checkMcf(mcf, mcf.run(),
   407              gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   408     checkMcf(mcf, mcf.lowerMap(l2).run(),
   409              gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   410     mcf.supplyMap(s5);
   411     checkMcf(mcf, mcf.run(),
   412              gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   413 
   414     // Check negative costs
   415     NetworkSimplex<Digraph> neg_mcf(neg_gr);
   416     neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
   417     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
   418       neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
   419     neg_mcf.upperMap(neg_u2);
   420     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
   421       neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
   422     neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
   423     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
   424       neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");
   425 
   426     NetworkSimplex<Digraph> negs_mcf(negs_gr);
   427     negs_mcf.costMap(negs_c).supplyMap(negs_s);
   428     checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
   429       negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);
   430   }
   431 
   432   // B. Test NetworkSimplex with each pivot rule
   433   {
   434     NetworkSimplex<Digraph> mcf(gr);
   435     mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
   436 
   437     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
   438              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
   439     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
   440              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
   441     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
   442              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
   443     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
   444              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
   445     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
   446              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
   447   }
   448 
   449   return 0;
   450 }