test/min_cost_flow_test.cc
author Alpar Juttner <alpar@cs.elte.hu>
Fri, 08 May 2009 16:21:06 +0100
changeset 655 c706534d4740
parent 640 6c408d864fa1
child 664 cc61d09f053b
permissions -rw-r--r--
Fix the faulty merge of unify-sources.sh in [586b65073025]
     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 )
   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     opt = (sum == supply[n]) || (pi[n] == 0);
   197   }
   198   
   199   return opt;
   200 }
   201 
   202 // Run a minimum cost flow algorithm and check the results
   203 template < typename MCF, typename GR,
   204            typename LM, typename UM,
   205            typename CM, typename SM,
   206            typename PT >
   207 void checkMcf( const MCF& mcf, PT mcf_result,
   208                const GR& gr, const LM& lower, const UM& upper,
   209                const CM& cost, const SM& supply,
   210                PT result, bool optimal, typename CM::Value total,
   211                const std::string &test_id = "",
   212                SupplyType type = EQ )
   213 {
   214   check(mcf_result == result, "Wrong result " + test_id);
   215   if (optimal) {
   216     typename GR::template ArcMap<typename SM::Value> flow(gr);
   217     typename GR::template NodeMap<typename CM::Value> pi(gr);
   218     mcf.flowMap(flow);
   219     mcf.potentialMap(pi);
   220     check(checkFlow(gr, lower, upper, supply, flow, type),
   221           "The flow is not feasible " + test_id);
   222     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   223     check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
   224           "Wrong potentials " + test_id);
   225   }
   226 }
   227 
   228 int main()
   229 {
   230   // Check the interfaces
   231   {
   232     typedef concepts::Digraph GR;
   233     checkConcept< McfClassConcept<GR, int, int>,
   234                   NetworkSimplex<GR> >();
   235     checkConcept< McfClassConcept<GR, double, double>,
   236                   NetworkSimplex<GR, double> >();
   237     checkConcept< McfClassConcept<GR, int, double>,
   238                   NetworkSimplex<GR, int, double> >();
   239   }
   240 
   241   // Run various MCF tests
   242   typedef ListDigraph Digraph;
   243   DIGRAPH_TYPEDEFS(ListDigraph);
   244 
   245   // Read the test digraph
   246   Digraph gr;
   247   Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
   248   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
   249   ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   250   Node v, w;
   251 
   252   std::istringstream input(test_lgf);
   253   DigraphReader<Digraph>(gr, input)
   254     .arcMap("cost", c)
   255     .arcMap("cap", u)
   256     .arcMap("low1", l1)
   257     .arcMap("low2", l2)
   258     .arcMap("low3", l3)
   259     .nodeMap("sup1", s1)
   260     .nodeMap("sup2", s2)
   261     .nodeMap("sup3", s3)
   262     .nodeMap("sup4", s4)
   263     .nodeMap("sup5", s5)
   264     .nodeMap("sup6", s6)
   265     .node("source", v)
   266     .node("target", w)
   267     .run();
   268   
   269   // Build a test digraph for testing negative costs
   270   Digraph ngr;
   271   Node n1 = ngr.addNode();
   272   Node n2 = ngr.addNode();
   273   Node n3 = ngr.addNode();
   274   Node n4 = ngr.addNode();
   275   Node n5 = ngr.addNode();
   276   Node n6 = ngr.addNode();
   277   Node n7 = ngr.addNode();
   278   
   279   Arc a1 = ngr.addArc(n1, n2);
   280   Arc a2 = ngr.addArc(n1, n3);
   281   Arc a3 = ngr.addArc(n2, n4);
   282   Arc a4 = ngr.addArc(n3, n4);
   283   Arc a5 = ngr.addArc(n3, n2);
   284   Arc a6 = ngr.addArc(n5, n3);
   285   Arc a7 = ngr.addArc(n5, n6);
   286   Arc a8 = ngr.addArc(n6, n7);
   287   Arc a9 = ngr.addArc(n7, n5);
   288   
   289   Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
   290   ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
   291   Digraph::NodeMap<int> ns(ngr, 0);
   292   
   293   nl2[a7] =  1000;
   294   nl2[a8] = -1000;
   295   
   296   ns[n1] =  100;
   297   ns[n4] = -100;
   298   
   299   nc[a1] =  100;
   300   nc[a2] =   30;
   301   nc[a3] =   20;
   302   nc[a4] =   80;
   303   nc[a5] =   50;
   304   nc[a6] =   10;
   305   nc[a7] =   80;
   306   nc[a8] =   30;
   307   nc[a9] = -120;
   308 
   309   // A. Test NetworkSimplex with the default pivot rule
   310   {
   311     NetworkSimplex<Digraph> mcf(gr);
   312 
   313     // Check the equality form
   314     mcf.upperMap(u).costMap(c);
   315     checkMcf(mcf, mcf.supplyMap(s1).run(),
   316              gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
   317     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   318              gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
   319     mcf.lowerMap(l2);
   320     checkMcf(mcf, mcf.supplyMap(s1).run(),
   321              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
   322     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   323              gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
   324     mcf.reset();
   325     checkMcf(mcf, mcf.supplyMap(s1).run(),
   326              gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
   327     checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
   328              gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
   329     mcf.reset();
   330     checkMcf(mcf, mcf.run(),
   331              gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
   332     checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
   333              gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
   334     mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
   335     checkMcf(mcf, mcf.run(),
   336              gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
   337 
   338     // Check the GEQ form
   339     mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
   340     checkMcf(mcf, mcf.run(),
   341              gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
   342     mcf.supplyType(mcf.GEQ);
   343     checkMcf(mcf, mcf.lowerMap(l2).run(),
   344              gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   345     mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
   346     checkMcf(mcf, mcf.run(),
   347              gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   348 
   349     // Check the LEQ form
   350     mcf.reset().supplyType(mcf.LEQ);
   351     mcf.upperMap(u).costMap(c).supplyMap(s6);
   352     checkMcf(mcf, mcf.run(),
   353              gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   354     checkMcf(mcf, mcf.lowerMap(l2).run(),
   355              gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   356     mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
   357     checkMcf(mcf, mcf.run(),
   358              gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   359 
   360     // Check negative costs
   361     NetworkSimplex<Digraph> nmcf(ngr);
   362     nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
   363     checkMcf(nmcf, nmcf.run(),
   364       ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
   365     checkMcf(nmcf, nmcf.upperMap(nu2).run(),
   366       ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
   367     nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
   368     checkMcf(nmcf, nmcf.run(),
   369       ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
   370   }
   371 
   372   // B. Test NetworkSimplex with each pivot rule
   373   {
   374     NetworkSimplex<Digraph> mcf(gr);
   375     mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
   376 
   377     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
   378              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
   379     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
   380              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
   381     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
   382              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
   383     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
   384              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
   385     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
   386              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
   387   }
   388 
   389   return 0;
   390 }