test/min_cost_flow_test.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 24 Apr 2009 12:22:06 +0200
changeset 613 b1811c363299
parent 607 9ad8d2122b50
child 615 e3d9bff447ed
permissions -rw-r--r--
Bug fix in NetworkSimplex (#234)
     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 
    22 #include <lemon/list_graph.h>
    23 #include <lemon/lgf_reader.h>
    24 
    25 #include <lemon/network_simplex.h>
    26 
    27 #include <lemon/concepts/digraph.h>
    28 #include <lemon/concept_check.h>
    29 
    30 #include "test_tools.h"
    31 
    32 using namespace lemon;
    33 
    34 char test_lgf[] =
    35   "@nodes\n"
    36   "label  sup1 sup2 sup3 sup4 sup5\n"
    37   "    1    20   27    0   20   30\n"
    38   "    2    -4    0    0   -8   -3\n"
    39   "    3     0    0    0    0    0\n"
    40   "    4     0    0    0    0    0\n"
    41   "    5     9    0    0    6   11\n"
    42   "    6    -6    0    0   -5   -6\n"
    43   "    7     0    0    0    0    0\n"
    44   "    8     0    0    0    0    3\n"
    45   "    9     3    0    0    0    0\n"
    46   "   10    -2    0    0   -7   -2\n"
    47   "   11     0    0    0  -10    0\n"
    48   "   12   -20  -27    0  -30  -20\n"
    49   "\n"
    50   "@arcs\n"
    51   "       cost  cap low1 low2\n"
    52   " 1  2    70   11    0    8\n"
    53   " 1  3   150    3    0    1\n"
    54   " 1  4    80   15    0    2\n"
    55   " 2  8    80   12    0    0\n"
    56   " 3  5   140    5    0    3\n"
    57   " 4  6    60   10    0    1\n"
    58   " 4  7    80    2    0    0\n"
    59   " 4  8   110    3    0    0\n"
    60   " 5  7    60   14    0    0\n"
    61   " 5 11   120   12    0    0\n"
    62   " 6  3     0    3    0    0\n"
    63   " 6  9   140    4    0    0\n"
    64   " 6 10    90    8    0    0\n"
    65   " 7  1    30    5    0    0\n"
    66   " 8 12    60   16    0    4\n"
    67   " 9 12    50    6    0    0\n"
    68   "10 12    70   13    0    5\n"
    69   "10  2   100    7    0    0\n"
    70   "10  7    60   10    0    0\n"
    71   "11 10    20   14    0    6\n"
    72   "12 11    30   10    0    0\n"
    73   "\n"
    74   "@attributes\n"
    75   "source 1\n"
    76   "target 12\n";
    77 
    78 
    79 enum ProblemType {
    80   EQ,
    81   GEQ,
    82   LEQ
    83 };
    84 
    85 // Check the interface of an MCF algorithm
    86 template <typename GR, typename Flow, typename Cost>
    87 class McfClassConcept
    88 {
    89 public:
    90 
    91   template <typename MCF>
    92   struct Constraints {
    93     void constraints() {
    94       checkConcept<concepts::Digraph, GR>();
    95 
    96       MCF mcf(g);
    97 
    98       b = mcf.reset()
    99              .lowerMap(lower)
   100              .upperMap(upper)
   101              .capacityMap(upper)
   102              .boundMaps(lower, upper)
   103              .costMap(cost)
   104              .supplyMap(sup)
   105              .stSupply(n, n, k)
   106              .flowMap(flow)
   107              .potentialMap(pot)
   108              .run();
   109       
   110       const MCF& const_mcf = mcf;
   111 
   112       const typename MCF::FlowMap &fm = const_mcf.flowMap();
   113       const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
   114 
   115       v = const_mcf.totalCost();
   116       double x = const_mcf.template totalCost<double>();
   117       v = const_mcf.flow(a);
   118       v = const_mcf.potential(n);
   119 
   120       ignore_unused_variable_warning(fm);
   121       ignore_unused_variable_warning(pm);
   122       ignore_unused_variable_warning(x);
   123     }
   124 
   125     typedef typename GR::Node Node;
   126     typedef typename GR::Arc Arc;
   127     typedef concepts::ReadMap<Node, Flow> NM;
   128     typedef concepts::ReadMap<Arc, Flow> FAM;
   129     typedef concepts::ReadMap<Arc, Cost> CAM;
   130 
   131     const GR &g;
   132     const FAM &lower;
   133     const FAM &upper;
   134     const CAM &cost;
   135     const NM &sup;
   136     const Node &n;
   137     const Arc &a;
   138     const Flow &k;
   139     Flow v;
   140     bool b;
   141 
   142     typename MCF::FlowMap &flow;
   143     typename MCF::PotentialMap &pot;
   144   };
   145 
   146 };
   147 
   148 
   149 // Check the feasibility of the given flow (primal soluiton)
   150 template < typename GR, typename LM, typename UM,
   151            typename SM, typename FM >
   152 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
   153                 const SM& supply, const FM& flow,
   154                 ProblemType type = EQ )
   155 {
   156   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   157 
   158   for (ArcIt e(gr); e != INVALID; ++e) {
   159     if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
   160   }
   161 
   162   for (NodeIt n(gr); n != INVALID; ++n) {
   163     typename SM::Value sum = 0;
   164     for (OutArcIt e(gr, n); e != INVALID; ++e)
   165       sum += flow[e];
   166     for (InArcIt e(gr, n); e != INVALID; ++e)
   167       sum -= flow[e];
   168     bool b = (type ==  EQ && sum == supply[n]) ||
   169              (type == GEQ && sum >= supply[n]) ||
   170              (type == LEQ && sum <= supply[n]);
   171     if (!b) return false;
   172   }
   173 
   174   return true;
   175 }
   176 
   177 // Check the feasibility of the given potentials (dual soluiton)
   178 // using the "Complementary Slackness" optimality condition
   179 template < typename GR, typename LM, typename UM,
   180            typename CM, typename SM, typename FM, typename PM >
   181 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   182                      const CM& cost, const SM& supply, const FM& flow, 
   183                      const PM& pi )
   184 {
   185   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   186 
   187   bool opt = true;
   188   for (ArcIt e(gr); opt && e != INVALID; ++e) {
   189     typename CM::Value red_cost =
   190       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
   191     opt = red_cost == 0 ||
   192           (red_cost > 0 && flow[e] == lower[e]) ||
   193           (red_cost < 0 && flow[e] == upper[e]);
   194   }
   195   
   196   for (NodeIt n(gr); opt && n != INVALID; ++n) {
   197     typename SM::Value sum = 0;
   198     for (OutArcIt e(gr, n); e != INVALID; ++e)
   199       sum += flow[e];
   200     for (InArcIt e(gr, n); e != INVALID; ++e)
   201       sum -= flow[e];
   202     opt = (sum == supply[n]) || (pi[n] == 0);
   203   }
   204   
   205   return opt;
   206 }
   207 
   208 // Run a minimum cost flow algorithm and check the results
   209 template < typename MCF, typename GR,
   210            typename LM, typename UM,
   211            typename CM, typename SM >
   212 void checkMcf( const MCF& mcf, bool mcf_result,
   213                const GR& gr, const LM& lower, const UM& upper,
   214                const CM& cost, const SM& supply,
   215                bool result, typename CM::Value total,
   216                const std::string &test_id = "",
   217                ProblemType type = EQ )
   218 {
   219   check(mcf_result == result, "Wrong result " + test_id);
   220   if (result) {
   221     check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
   222           "The flow is not feasible " + test_id);
   223     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   224     check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
   225                          mcf.potentialMap()),
   226           "Wrong potentials " + test_id);
   227   }
   228 }
   229 
   230 int main()
   231 {
   232   // Check the interfaces
   233   {
   234     typedef int Flow;
   235     typedef int Cost;
   236     // TODO: This typedef should be enabled if the standard maps are
   237     // reference maps in the graph concepts (See #190).
   238 /**/
   239     //typedef concepts::Digraph GR;
   240     typedef ListDigraph GR;
   241 /**/
   242     checkConcept< McfClassConcept<GR, Flow, Cost>,
   243                   NetworkSimplex<GR, Flow, Cost> >();
   244   }
   245 
   246   // Run various MCF tests
   247   typedef ListDigraph Digraph;
   248   DIGRAPH_TYPEDEFS(ListDigraph);
   249 
   250   // Read the test digraph
   251   Digraph gr;
   252   Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
   253   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
   254   ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   255   Node v, w;
   256 
   257   std::istringstream input(test_lgf);
   258   DigraphReader<Digraph>(gr, input)
   259     .arcMap("cost", c)
   260     .arcMap("cap", u)
   261     .arcMap("low1", l1)
   262     .arcMap("low2", l2)
   263     .nodeMap("sup1", s1)
   264     .nodeMap("sup2", s2)
   265     .nodeMap("sup3", s3)
   266     .nodeMap("sup4", s4)
   267     .nodeMap("sup5", s5)
   268     .node("source", v)
   269     .node("target", w)
   270     .run();
   271 
   272   // A. Test NetworkSimplex with the default pivot rule
   273   {
   274     NetworkSimplex<Digraph> mcf(gr);
   275 
   276     // Check the equality form
   277     mcf.upperMap(u).costMap(c);
   278     checkMcf(mcf, mcf.supplyMap(s1).run(),
   279              gr, l1, u, c, s1, true,  5240, "#A1");
   280     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   281              gr, l1, u, c, s2, true,  7620, "#A2");
   282     mcf.lowerMap(l2);
   283     checkMcf(mcf, mcf.supplyMap(s1).run(),
   284              gr, l2, u, c, s1, true,  5970, "#A3");
   285     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   286              gr, l2, u, c, s2, true,  8010, "#A4");
   287     mcf.reset();
   288     checkMcf(mcf, mcf.supplyMap(s1).run(),
   289              gr, l1, cu, cc, s1, true,  74, "#A5");
   290     checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
   291              gr, l2, cu, cc, s2, true,  94, "#A6");
   292     mcf.reset();
   293     checkMcf(mcf, mcf.run(),
   294              gr, l1, cu, cc, s3, true,   0, "#A7");
   295     checkMcf(mcf, mcf.boundMaps(l2, u).run(),
   296              gr, l2, u, cc, s3, false,   0, "#A8");
   297 
   298     // Check the GEQ form
   299     mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
   300     checkMcf(mcf, mcf.run(),
   301              gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
   302     mcf.problemType(mcf.GEQ);
   303     checkMcf(mcf, mcf.lowerMap(l2).run(),
   304              gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
   305     mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
   306     checkMcf(mcf, mcf.run(),
   307              gr, l2, u, c, s5, false,    0, "#A11", GEQ);
   308 
   309     // Check the LEQ form
   310     mcf.reset().problemType(mcf.LEQ);
   311     mcf.upperMap(u).costMap(c).supplyMap(s5);
   312     checkMcf(mcf, mcf.run(),
   313              gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
   314     checkMcf(mcf, mcf.lowerMap(l2).run(),
   315              gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
   316     mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
   317     checkMcf(mcf, mcf.run(),
   318              gr, l2, u, c, s4, false,    0, "#A14", LEQ);
   319   }
   320 
   321   // B. Test NetworkSimplex with each pivot rule
   322   {
   323     NetworkSimplex<Digraph> mcf(gr);
   324     mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
   325 
   326     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
   327              gr, l2, u, c, s1, true,  5970, "#B1");
   328     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
   329              gr, l2, u, c, s1, true,  5970, "#B2");
   330     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
   331              gr, l2, u, c, s1, true,  5970, "#B3");
   332     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
   333              gr, l2, u, c, s1, true,  5970, "#B4");
   334     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
   335              gr, l2, u, c, s1, true,  5970, "#B5");
   336   }
   337 
   338   return 0;
   339 }