test/min_cost_flow_test.cc
author Akos Ladanyi <ladanyi@tmit.bme.hu>
Tue, 28 Apr 2009 16:49:57 +0100
changeset 630 936355864d6e
parent 609 e6927fe719e6
child 640 6c408d864fa1
permissions -rw-r--r--
Add cbc.h to lemon/Makefile.am (#279)
     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     typedef concepts::Digraph GR;
   237     checkConcept< McfClassConcept<GR, Flow, Cost>,
   238                   NetworkSimplex<GR, Flow, Cost> >();
   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), u(gr);
   248   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(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     .nodeMap("sup1", s1)
   259     .nodeMap("sup2", s2)
   260     .nodeMap("sup3", s3)
   261     .nodeMap("sup4", s4)
   262     .nodeMap("sup5", s5)
   263     .node("source", v)
   264     .node("target", w)
   265     .run();
   266 
   267   // A. Test NetworkSimplex with the default pivot rule
   268   {
   269     NetworkSimplex<Digraph> mcf(gr);
   270 
   271     // Check the equality form
   272     mcf.upperMap(u).costMap(c);
   273     checkMcf(mcf, mcf.supplyMap(s1).run(),
   274              gr, l1, u, c, s1, true,  5240, "#A1");
   275     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   276              gr, l1, u, c, s2, true,  7620, "#A2");
   277     mcf.lowerMap(l2);
   278     checkMcf(mcf, mcf.supplyMap(s1).run(),
   279              gr, l2, u, c, s1, true,  5970, "#A3");
   280     checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   281              gr, l2, u, c, s2, true,  8010, "#A4");
   282     mcf.reset();
   283     checkMcf(mcf, mcf.supplyMap(s1).run(),
   284              gr, l1, cu, cc, s1, true,  74, "#A5");
   285     checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
   286              gr, l2, cu, cc, s2, true,  94, "#A6");
   287     mcf.reset();
   288     checkMcf(mcf, mcf.run(),
   289              gr, l1, cu, cc, s3, true,   0, "#A7");
   290     checkMcf(mcf, mcf.boundMaps(l2, u).run(),
   291              gr, l2, u, cc, s3, false,   0, "#A8");
   292 
   293     // Check the GEQ form
   294     mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
   295     checkMcf(mcf, mcf.run(),
   296              gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
   297     mcf.problemType(mcf.GEQ);
   298     checkMcf(mcf, mcf.lowerMap(l2).run(),
   299              gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
   300     mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
   301     checkMcf(mcf, mcf.run(),
   302              gr, l2, u, c, s5, false,    0, "#A11", GEQ);
   303 
   304     // Check the LEQ form
   305     mcf.reset().problemType(mcf.LEQ);
   306     mcf.upperMap(u).costMap(c).supplyMap(s5);
   307     checkMcf(mcf, mcf.run(),
   308              gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
   309     checkMcf(mcf, mcf.lowerMap(l2).run(),
   310              gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
   311     mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
   312     checkMcf(mcf, mcf.run(),
   313              gr, l2, u, c, s4, false,    0, "#A14", LEQ);
   314   }
   315 
   316   // B. Test NetworkSimplex with each pivot rule
   317   {
   318     NetworkSimplex<Digraph> mcf(gr);
   319     mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
   320 
   321     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
   322              gr, l2, u, c, s1, true,  5970, "#B1");
   323     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
   324              gr, l2, u, c, s1, true,  5970, "#B2");
   325     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
   326              gr, l2, u, c, s1, true,  5970, "#B3");
   327     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
   328              gr, l2, u, c, s1, true,  5970, "#B4");
   329     checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
   330              gr, l2, u, c, s1, true,  5970, "#B5");
   331   }
   332 
   333   return 0;
   334 }