test/min_cost_flow_test.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Wed, 29 Apr 2009 03:15:24 +0200
changeset 687 6c408d864fa1
parent 662 e3d9bff447ed
child 689 111698359429
permissions -rw-r--r--
Support negative costs and bounds in NetworkSimplex (#270)

* The interface is reworked to support negative costs and bounds.
- ProblemType and problemType() are renamed to
SupplyType and supplyType(), see also #234.
- ProblemType type is introduced similarly to the LP interface.
- 'bool run()' is replaced by 'ProblemType run()' to handle
unbounded problem instances, as well.
- Add INF public member constant similarly to the LP interface.
* Remove capacityMap() and boundMaps(), see also #266.
* Update the problem definition in the MCF module.
* Remove the usage of Circulation (and adaptors) for checking feasibility.
Check feasibility by examining the artifical arcs instead (after solving
the problem).
* Additional check for unbounded negative cycles found during the
algorithm (it is possible now, since negative costs are allowed).
* Fix in the constructor (the value types needn't be integer any more),
see also #254.
* Improve and extend the doc.
* Rework the test file and add test cases for negative costs and bounds.
     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 Flow, 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 
    99       b = mcf.reset()
   100              .lowerMap(lower)
   101              .upperMap(upper)
   102              .costMap(cost)
   103              .supplyMap(sup)
   104              .stSupply(n, n, k)
   105              .flowMap(flow)
   106              .potentialMap(pot)
   107              .run();
   108       
   109       const MCF& const_mcf = mcf;
   110 
   111       const typename MCF::FlowMap &fm = const_mcf.flowMap();
   112       const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
   113 
   114       c = const_mcf.totalCost();
   115       double x = const_mcf.template totalCost<double>();
   116       v = const_mcf.flow(a);
   117       c = const_mcf.potential(n);
   118       
   119       v = const_mcf.INF;
   120 
   121       ignore_unused_variable_warning(fm);
   122       ignore_unused_variable_warning(pm);
   123       ignore_unused_variable_warning(x);
   124     }
   125 
   126     typedef typename GR::Node Node;
   127     typedef typename GR::Arc Arc;
   128     typedef concepts::ReadMap<Node, Flow> NM;
   129     typedef concepts::ReadMap<Arc, Flow> FAM;
   130     typedef concepts::ReadMap<Arc, Cost> CAM;
   131 
   132     const GR &g;
   133     const FAM &lower;
   134     const FAM &upper;
   135     const CAM &cost;
   136     const NM &sup;
   137     const Node &n;
   138     const Arc &a;
   139     const Flow &k;
   140     Flow v;
   141     Cost c;
   142     bool b;
   143 
   144     typename MCF::FlowMap &flow;
   145     typename MCF::PotentialMap &pot;
   146   };
   147 
   148 };
   149 
   150 
   151 // Check the feasibility of the given flow (primal soluiton)
   152 template < typename GR, typename LM, typename UM,
   153            typename SM, typename FM >
   154 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
   155                 const SM& supply, const FM& flow,
   156                 SupplyType type = EQ )
   157 {
   158   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   159 
   160   for (ArcIt e(gr); e != INVALID; ++e) {
   161     if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
   162   }
   163 
   164   for (NodeIt n(gr); n != INVALID; ++n) {
   165     typename SM::Value sum = 0;
   166     for (OutArcIt e(gr, n); e != INVALID; ++e)
   167       sum += flow[e];
   168     for (InArcIt e(gr, n); e != INVALID; ++e)
   169       sum -= flow[e];
   170     bool b = (type ==  EQ && sum == supply[n]) ||
   171              (type == GEQ && sum >= supply[n]) ||
   172              (type == LEQ && sum <= supply[n]);
   173     if (!b) return false;
   174   }
   175 
   176   return true;
   177 }
   178 
   179 // Check the feasibility of the given potentials (dual soluiton)
   180 // using the "Complementary Slackness" optimality condition
   181 template < typename GR, typename LM, typename UM,
   182            typename CM, typename SM, typename FM, typename PM >
   183 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   184                      const CM& cost, const SM& supply, const FM& flow, 
   185                      const PM& pi )
   186 {
   187   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   188 
   189   bool opt = true;
   190   for (ArcIt e(gr); opt && e != INVALID; ++e) {
   191     typename CM::Value red_cost =
   192       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
   193     opt = red_cost == 0 ||
   194           (red_cost > 0 && flow[e] == lower[e]) ||
   195           (red_cost < 0 && flow[e] == upper[e]);
   196   }
   197   
   198   for (NodeIt n(gr); opt && n != INVALID; ++n) {
   199     typename SM::Value sum = 0;
   200     for (OutArcIt e(gr, n); e != INVALID; ++e)
   201       sum += flow[e];
   202     for (InArcIt e(gr, n); e != INVALID; ++e)
   203       sum -= flow[e];
   204     opt = (sum == supply[n]) || (pi[n] == 0);
   205   }
   206   
   207   return opt;
   208 }
   209 
   210 // Run a minimum cost flow algorithm and check the results
   211 template < typename MCF, typename GR,
   212            typename LM, typename UM,
   213            typename CM, typename SM,
   214            typename PT >
   215 void checkMcf( const MCF& mcf, PT mcf_result,
   216                const GR& gr, const LM& lower, const UM& upper,
   217                const CM& cost, const SM& supply,
   218                PT result, bool optimal, typename CM::Value total,
   219                const std::string &test_id = "",
   220                SupplyType type = EQ )
   221 {
   222   check(mcf_result == result, "Wrong result " + test_id);
   223   if (optimal) {
   224     check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
   225           "The flow is not feasible " + test_id);
   226     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   227     check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
   228                          mcf.potentialMap()),
   229           "Wrong potentials " + test_id);
   230   }
   231 }
   232 
   233 int main()
   234 {
   235   // Check the interfaces
   236   {
   237     typedef int Flow;
   238     typedef int Cost;
   239     typedef concepts::Digraph GR;
   240     checkConcept< McfClassConcept<GR, Flow, Cost>,
   241                   NetworkSimplex<GR, Flow, Cost> >();
   242   }
   243 
   244   // Run various MCF tests
   245   typedef ListDigraph Digraph;
   246   DIGRAPH_TYPEDEFS(ListDigraph);
   247 
   248   // Read the test digraph
   249   Digraph gr;
   250   Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
   251   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
   252   ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   253   Node v, w;
   254 
   255   std::istringstream input(test_lgf);
   256   DigraphReader<Digraph>(gr, input)
   257     .arcMap("cost", c)
   258     .arcMap("cap", u)
   259     .arcMap("low1", l1)
   260     .arcMap("low2", l2)
   261     .arcMap("low3", l3)
   262     .nodeMap("sup1", s1)
   263     .nodeMap("sup2", s2)
   264     .nodeMap("sup3", s3)
   265     .nodeMap("sup4", s4)
   266     .nodeMap("sup5", s5)
   267     .nodeMap("sup6", s6)
   268     .node("source", v)
   269     .node("target", w)
   270     .run();
   271   
   272   // Build a test digraph for testing negative costs
   273   Digraph ngr;
   274   Node n1 = ngr.addNode();
   275   Node n2 = ngr.addNode();
   276   Node n3 = ngr.addNode();
   277   Node n4 = ngr.addNode();
   278   Node n5 = ngr.addNode();
   279   Node n6 = ngr.addNode();
   280   Node n7 = ngr.addNode();
   281   
   282   Arc a1 = ngr.addArc(n1, n2);
   283   Arc a2 = ngr.addArc(n1, n3);
   284   Arc a3 = ngr.addArc(n2, n4);
   285   Arc a4 = ngr.addArc(n3, n4);
   286   Arc a5 = ngr.addArc(n3, n2);
   287   Arc a6 = ngr.addArc(n5, n3);
   288   Arc a7 = ngr.addArc(n5, n6);
   289   Arc a8 = ngr.addArc(n6, n7);
   290   Arc a9 = ngr.addArc(n7, n5);
   291   
   292   Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
   293   ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
   294   Digraph::NodeMap<int> ns(ngr, 0);
   295   
   296   nl2[a7] =  1000;
   297   nl2[a8] = -1000;
   298   
   299   ns[n1] =  100;
   300   ns[n4] = -100;
   301   
   302   nc[a1] =  100;
   303   nc[a2] =   30;
   304   nc[a3] =   20;
   305   nc[a4] =   80;
   306   nc[a5] =   50;
   307   nc[a6] =   10;
   308   nc[a7] =   80;
   309   nc[a8] =   30;
   310   nc[a9] = -120;
   311 
   312   // A. Test NetworkSimplex with the default pivot rule
   313   {
   314     NetworkSimplex<Digraph> mcf(gr);
   315 
   316     // Check the equality form
   317     mcf.upperMap(u).costMap(c);
   318     checkMcf(mcf, mcf.supplyMap(s1).run(),
   319              gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
   320     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   321              gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
   322     mcf.lowerMap(l2);
   323     checkMcf(mcf, mcf.supplyMap(s1).run(),
   324              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
   325     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   326              gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
   327     mcf.reset();
   328     checkMcf(mcf, mcf.supplyMap(s1).run(),
   329              gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
   330     checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
   331              gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
   332     mcf.reset();
   333     checkMcf(mcf, mcf.run(),
   334              gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
   335     checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
   336              gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
   337     mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
   338     checkMcf(mcf, mcf.run(),
   339              gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
   340 
   341     // Check the GEQ form
   342     mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
   343     checkMcf(mcf, mcf.run(),
   344              gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
   345     mcf.supplyType(mcf.GEQ);
   346     checkMcf(mcf, mcf.lowerMap(l2).run(),
   347              gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   348     mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
   349     checkMcf(mcf, mcf.run(),
   350              gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   351 
   352     // Check the LEQ form
   353     mcf.reset().supplyType(mcf.LEQ);
   354     mcf.upperMap(u).costMap(c).supplyMap(s6);
   355     checkMcf(mcf, mcf.run(),
   356              gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   357     checkMcf(mcf, mcf.lowerMap(l2).run(),
   358              gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   359     mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
   360     checkMcf(mcf, mcf.run(),
   361              gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   362 
   363     // Check negative costs
   364     NetworkSimplex<Digraph> nmcf(ngr);
   365     nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
   366     checkMcf(nmcf, nmcf.run(),
   367       ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
   368     checkMcf(nmcf, nmcf.upperMap(nu2).run(),
   369       ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
   370     nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
   371     checkMcf(nmcf, nmcf.run(),
   372       ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
   373   }
   374 
   375   // B. Test NetworkSimplex with each pivot rule
   376   {
   377     NetworkSimplex<Digraph> mcf(gr);
   378     mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
   379 
   380     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
   381              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
   382     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
   383              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
   384     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
   385              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
   386     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
   387              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
   388     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
   389              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
   390   }
   391 
   392   return 0;
   393 }