test/min_cost_flow_test.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Tue, 12 May 2009 12:08:06 +0200
changeset 711 cc61d09f053b
parent 689 111698359429
child 716 4faca85d40e6
permissions -rw-r--r--
Extend min cost flow test file + check dual costs (#291)
     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 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       MCF mcf(g);
    98       const MCF& const_mcf = mcf;
    99 
   100       b = mcf.reset()
   101              .lowerMap(lower)
   102              .upperMap(upper)
   103              .costMap(cost)
   104              .supplyMap(sup)
   105              .stSupply(n, n, k)
   106              .run();
   107 
   108       c = const_mcf.totalCost();
   109       x = const_mcf.template totalCost<double>();
   110       v = const_mcf.flow(a);
   111       c = const_mcf.potential(n);
   112       const_mcf.flowMap(fm);
   113       const_mcf.potentialMap(pm);
   114     }
   115 
   116     typedef typename GR::Node Node;
   117     typedef typename GR::Arc Arc;
   118     typedef concepts::ReadMap<Node, Value> NM;
   119     typedef concepts::ReadMap<Arc, Value> VAM;
   120     typedef concepts::ReadMap<Arc, Cost> CAM;
   121     typedef concepts::WriteMap<Arc, Value> FlowMap;
   122     typedef concepts::WriteMap<Node, Cost> PotMap;
   123 
   124     const GR &g;
   125     const VAM &lower;
   126     const VAM &upper;
   127     const CAM &cost;
   128     const NM &sup;
   129     const Node &n;
   130     const Arc &a;
   131     const Value &k;
   132     FlowMap fm;
   133     PotMap pm;
   134     bool b;
   135     double x;
   136     typename MCF::Value v;
   137     typename MCF::Cost c;
   138   };
   139 
   140 };
   141 
   142 
   143 // Check the feasibility of the given flow (primal soluiton)
   144 template < typename GR, typename LM, typename UM,
   145            typename SM, typename FM >
   146 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
   147                 const SM& supply, const FM& flow,
   148                 SupplyType type = EQ )
   149 {
   150   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   151 
   152   for (ArcIt e(gr); e != INVALID; ++e) {
   153     if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
   154   }
   155 
   156   for (NodeIt n(gr); n != INVALID; ++n) {
   157     typename SM::Value sum = 0;
   158     for (OutArcIt e(gr, n); e != INVALID; ++e)
   159       sum += flow[e];
   160     for (InArcIt e(gr, n); e != INVALID; ++e)
   161       sum -= flow[e];
   162     bool b = (type ==  EQ && sum == supply[n]) ||
   163              (type == GEQ && sum >= supply[n]) ||
   164              (type == LEQ && sum <= supply[n]);
   165     if (!b) return false;
   166   }
   167 
   168   return true;
   169 }
   170 
   171 // Check the feasibility of the given potentials (dual soluiton)
   172 // using the "Complementary Slackness" optimality condition
   173 template < typename GR, typename LM, typename UM,
   174            typename CM, typename SM, typename FM, typename PM >
   175 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   176                      const CM& cost, const SM& supply, const FM& flow, 
   177                      const PM& pi, SupplyType type )
   178 {
   179   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   180 
   181   bool opt = true;
   182   for (ArcIt e(gr); opt && e != INVALID; ++e) {
   183     typename CM::Value red_cost =
   184       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
   185     opt = red_cost == 0 ||
   186           (red_cost > 0 && flow[e] == lower[e]) ||
   187           (red_cost < 0 && flow[e] == upper[e]);
   188   }
   189   
   190   for (NodeIt n(gr); opt && n != INVALID; ++n) {
   191     typename SM::Value sum = 0;
   192     for (OutArcIt e(gr, n); e != INVALID; ++e)
   193       sum += flow[e];
   194     for (InArcIt e(gr, n); e != INVALID; ++e)
   195       sum -= flow[e];
   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     }
   201   }
   202   
   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;
   238 }
   239 
   240 // Run a minimum cost flow algorithm and check the results
   241 template < typename MCF, typename GR,
   242            typename LM, typename UM,
   243            typename CM, typename SM,
   244            typename PT >
   245 void checkMcf( const MCF& mcf, PT mcf_result,
   246                const GR& gr, const LM& lower, const UM& upper,
   247                const CM& cost, const SM& supply,
   248                PT result, bool optimal, typename CM::Value total,
   249                const std::string &test_id = "",
   250                SupplyType type = EQ )
   251 {
   252   check(mcf_result == result, "Wrong result " + test_id);
   253   if (optimal) {
   254     typename GR::template ArcMap<typename SM::Value> flow(gr);
   255     typename GR::template NodeMap<typename CM::Value> pi(gr);
   256     mcf.flowMap(flow);
   257     mcf.potentialMap(pi);
   258     check(checkFlow(gr, lower, upper, supply, flow, type),
   259           "The flow is not feasible " + test_id);
   260     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   261     check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
   262           "Wrong potentials " + test_id);
   263     check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
   264           "Wrong dual cost " + test_id);
   265   }
   266 }
   267 
   268 int main()
   269 {
   270   // Check the interfaces
   271   {
   272     typedef concepts::Digraph GR;
   273     checkConcept< McfClassConcept<GR, int, int>,
   274                   NetworkSimplex<GR> >();
   275     checkConcept< McfClassConcept<GR, double, double>,
   276                   NetworkSimplex<GR, double> >();
   277     checkConcept< McfClassConcept<GR, int, double>,
   278                   NetworkSimplex<GR, int, double> >();
   279   }
   280 
   281   // Run various MCF tests
   282   typedef ListDigraph Digraph;
   283   DIGRAPH_TYPEDEFS(ListDigraph);
   284 
   285   // Read the test digraph
   286   Digraph gr;
   287   Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
   288   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
   289   ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   290   Node v, w;
   291 
   292   std::istringstream input(test_lgf);
   293   DigraphReader<Digraph>(gr, input)
   294     .arcMap("cost", c)
   295     .arcMap("cap", u)
   296     .arcMap("low1", l1)
   297     .arcMap("low2", l2)
   298     .arcMap("low3", l3)
   299     .nodeMap("sup1", s1)
   300     .nodeMap("sup2", s2)
   301     .nodeMap("sup3", s3)
   302     .nodeMap("sup4", s4)
   303     .nodeMap("sup5", s5)
   304     .nodeMap("sup6", s6)
   305     .node("source", v)
   306     .node("target", w)
   307     .run();
   308   
   309   // Build test digraphs with negative costs
   310   Digraph neg_gr;
   311   Node n1 = neg_gr.addNode();
   312   Node n2 = neg_gr.addNode();
   313   Node n3 = neg_gr.addNode();
   314   Node n4 = neg_gr.addNode();
   315   Node n5 = neg_gr.addNode();
   316   Node n6 = neg_gr.addNode();
   317   Node n7 = neg_gr.addNode();
   318   
   319   Arc a1 = neg_gr.addArc(n1, n2);
   320   Arc a2 = neg_gr.addArc(n1, n3);
   321   Arc a3 = neg_gr.addArc(n2, n4);
   322   Arc a4 = neg_gr.addArc(n3, n4);
   323   Arc a5 = neg_gr.addArc(n3, n2);
   324   Arc a6 = neg_gr.addArc(n5, n3);
   325   Arc a7 = neg_gr.addArc(n5, n6);
   326   Arc a8 = neg_gr.addArc(n6, n7);
   327   Arc a9 = neg_gr.addArc(n7, n5);
   328   
   329   Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
   330   ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
   331   Digraph::NodeMap<int> neg_s(neg_gr, 0);
   332   
   333   neg_l2[a7] =  1000;
   334   neg_l2[a8] = -1000;
   335   
   336   neg_s[n1] =  100;
   337   neg_s[n4] = -100;
   338   
   339   neg_c[a1] =  100;
   340   neg_c[a2] =   30;
   341   neg_c[a3] =   20;
   342   neg_c[a4] =   80;
   343   neg_c[a5] =   50;
   344   neg_c[a6] =   10;
   345   neg_c[a7] =   80;
   346   neg_c[a8] =   30;
   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 
   359 
   360   // A. Test NetworkSimplex with the default pivot rule
   361   {
   362     NetworkSimplex<Digraph> mcf(gr);
   363 
   364     // Check the equality form
   365     mcf.upperMap(u).costMap(c);
   366     checkMcf(mcf, mcf.supplyMap(s1).run(),
   367              gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
   368     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   369              gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
   370     mcf.lowerMap(l2);
   371     checkMcf(mcf, mcf.supplyMap(s1).run(),
   372              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
   373     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   374              gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
   375     mcf.reset();
   376     checkMcf(mcf, mcf.supplyMap(s1).run(),
   377              gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
   378     checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
   379              gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
   380     mcf.reset();
   381     checkMcf(mcf, mcf.run(),
   382              gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
   383     checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
   384              gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
   385     mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
   386     checkMcf(mcf, mcf.run(),
   387              gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
   388 
   389     // Check the GEQ form
   390     mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
   391     checkMcf(mcf, mcf.run(),
   392              gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
   393     mcf.supplyType(mcf.GEQ);
   394     checkMcf(mcf, mcf.lowerMap(l2).run(),
   395              gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   396     mcf.supplyMap(s6);
   397     checkMcf(mcf, mcf.run(),
   398              gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   399 
   400     // Check the LEQ form
   401     mcf.reset().supplyType(mcf.LEQ);
   402     mcf.upperMap(u).costMap(c).supplyMap(s6);
   403     checkMcf(mcf, mcf.run(),
   404              gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   405     checkMcf(mcf, mcf.lowerMap(l2).run(),
   406              gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   407     mcf.supplyMap(s5);
   408     checkMcf(mcf, mcf.run(),
   409              gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   410 
   411     // Check negative costs
   412     NetworkSimplex<Digraph> neg_mcf(neg_gr);
   413     neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
   414     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
   415       neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
   416     neg_mcf.upperMap(neg_u2);
   417     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
   418       neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
   419     neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
   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);
   427   }
   428 
   429   // B. Test NetworkSimplex with each pivot rule
   430   {
   431     NetworkSimplex<Digraph> mcf(gr);
   432     mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
   433 
   434     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
   435              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
   436     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
   437              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
   438     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
   439              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
   440     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
   441              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
   442     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
   443              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
   444   }
   445 
   446   return 0;
   447 }