lemon: test/min_cost_flow_test.cc@eb2f9d453070

     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-2010

     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  */

    19 #include <iostream>

    20 #include <fstream>

    21 #include <limits>

    23 #include <lemon/list_graph.h>

    24 #include <lemon/lgf_reader.h>

    26 #include <lemon/network_simplex.h>

    27 #include <lemon/capacity_scaling.h>

    28 #include <lemon/cost_scaling.h>

    29 #include <lemon/cycle_canceling.h>

    31 #include <lemon/concepts/digraph.h>

    32 #include <lemon/concepts/heap.h>

    33 #include <lemon/concept_check.h>

    35 #include "test_tools.h"

    37 using namespace lemon;

    39 // Test networks

    40 char test_lgf[] =

    41   "@nodes\n"

    42   "label  sup1 sup2 sup3 sup4 sup5 sup6\n"

    43   "    1    20   27    0   30   20   30\n"

    44   "    2    -4    0    0    0   -8   -3\n"

    45   "    3     0    0    0    0    0    0\n"

    46   "    4     0    0    0    0    0    0\n"

    47   "    5     9    0    0    0    6   11\n"

    48   "    6    -6    0    0    0   -5   -6\n"

    49   "    7     0    0    0    0    0    0\n"

    50   "    8     0    0    0    0    0    3\n"

    51   "    9     3    0    0    0    0    0\n"

    52   "   10    -2    0    0    0   -7   -2\n"

    53   "   11     0    0    0    0  -10    0\n"

    54   "   12   -20  -27    0  -30  -30  -20\n"

    55   "\n"

    56   "@arcs\n"

    57   "       cost  cap low1 low2 low3\n"

    58   " 1  2    70   11    0    8    8\n"

    59   " 1  3   150    3    0    1    0\n"

    60   " 1  4    80   15    0    2    2\n"

    61   " 2  8    80   12    0    0    0\n"

    62   " 3  5   140    5    0    3    1\n"

    63   " 4  6    60   10    0    1    0\n"

    64   " 4  7    80    2    0    0    0\n"

    65   " 4  8   110    3    0    0    0\n"

    66   " 5  7    60   14    0    0    0\n"

    67   " 5 11   120   12    0    0    0\n"

    68   " 6  3     0    3    0    0    0\n"

    69   " 6  9   140    4    0    0    0\n"

    70   " 6 10    90    8    0    0    0\n"

    71   " 7  1    30    5    0    0   -5\n"

    72   " 8 12    60   16    0    4    3\n"

    73   " 9 12    50    6    0    0    0\n"

    74   "10 12    70   13    0    5    2\n"

    75   "10  2   100    7    0    0    0\n"

    76   "10  7    60   10    0    0   -3\n"

    77   "11 10    20   14    0    6  -20\n"

    78   "12 11    30   10    0    0  -10\n"

    79   "\n"

    80   "@attributes\n"

    81   "source 1\n"

    82   "target 12\n";

    84 char test_neg1_lgf[] =

    85   "@nodes\n"

    86   "label   sup\n"

    87   "    1   100\n"

    88   "    2     0\n"

    89   "    3     0\n"

    90   "    4  -100\n"

    91   "    5     0\n"

    92   "    6     0\n"

    93   "    7     0\n"

    94   "@arcs\n"

    95   "      cost   low1   low2\n"

    96   "1 2    100      0      0\n"

    97   "1 3     30      0      0\n"

    98   "2 4     20      0      0\n"

    99   "3 4     80      0      0\n"

   100   "3 2     50      0      0\n"

   101   "5 3     10      0      0\n"

   102   "5 6     80      0   1000\n"

   103   "6 7     30      0  -1000\n"

   104   "7 5   -120      0      0\n";

   106 char test_neg2_lgf[] =

   107   "@nodes\n"

   108   "label   sup\n"

   109   "    1   100\n"

   110   "    2  -300\n"

   111   "@arcs\n"

   112   "      cost\n"

   113   "1 2     -1\n";

   116 // Test data

   117 typedef ListDigraph Digraph;

   118 DIGRAPH_TYPEDEFS(ListDigraph);

   120 Digraph gr;

   121 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);

   122 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);

   123 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());

   124 Node v, w;

   126 Digraph neg1_gr;

   127 Digraph::ArcMap<int> neg1_c(neg1_gr), neg1_l1(neg1_gr), neg1_l2(neg1_gr);

   128 ConstMap<Arc, int> neg1_u1(std::numeric_limits<int>::max()), neg1_u2(5000);

   129 Digraph::NodeMap<int> neg1_s(neg1_gr);

   131 Digraph neg2_gr;

   132 Digraph::ArcMap<int> neg2_c(neg2_gr);

   133 ConstMap<Arc, int> neg2_l(0), neg2_u(1000);

   134 Digraph::NodeMap<int> neg2_s(neg2_gr);

   137 enum SupplyType {

   138   EQ,

   139   GEQ,

   140   LEQ

   141 };

   144 // Check the interface of an MCF algorithm

   145 template <typename GR, typename Value, typename Cost>

   146 class McfClassConcept

   147 {

   148 public:

   150   template <typename MCF>

   151   struct Constraints {

   152     void constraints() {

   153       checkConcept<concepts::Digraph, GR>();

   155       const Constraints& me = *this;

   157       MCF mcf(me.g);

   158       const MCF& const_mcf = mcf;

   160       b = mcf.reset().resetParams()

   161              .lowerMap(me.lower)

   162              .upperMap(me.upper)

   163              .costMap(me.cost)

   164              .supplyMap(me.sup)

   165              .stSupply(me.n, me.n, me.k)

   166              .run();

   168       c = const_mcf.totalCost();

   169       x = const_mcf.template totalCost<double>();

   170       v = const_mcf.flow(me.a);

   171       c = const_mcf.potential(me.n);

   172       const_mcf.flowMap(fm);

   173       const_mcf.potentialMap(pm);

   174     }

   176     typedef typename GR::Node Node;

   177     typedef typename GR::Arc Arc;

   178     typedef concepts::ReadMap<Node, Value> NM;

   179     typedef concepts::ReadMap<Arc, Value> VAM;

   180     typedef concepts::ReadMap<Arc, Cost> CAM;

   181     typedef concepts::WriteMap<Arc, Value> FlowMap;

   182     typedef concepts::WriteMap<Node, Cost> PotMap;

   184     GR g;

   185     VAM lower;

   186     VAM upper;

   187     CAM cost;

   188     NM sup;

   189     Node n;

   190     Arc a;

   191     Value k;

   193     FlowMap fm;

   194     PotMap pm;

   195     bool b;

   196     double x;

   197     typename MCF::Value v;

   198     typename MCF::Cost c;

   199   };

   201 };

   204 // Check the feasibility of the given flow (primal soluiton)

   205 template < typename GR, typename LM, typename UM,

   206            typename SM, typename FM >

   207 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,

   208                 const SM& supply, const FM& flow,

   209                 SupplyType type = EQ )

   210 {

   211   TEMPLATE_DIGRAPH_TYPEDEFS(GR);

   213   for (ArcIt e(gr); e != INVALID; ++e) {

   214     if (flow[e] < lower[e] || flow[e] > upper[e]) return false;

   215   }

   217   for (NodeIt n(gr); n != INVALID; ++n) {

   218     typename SM::Value sum = 0;

   219     for (OutArcIt e(gr, n); e != INVALID; ++e)

   220       sum += flow[e];

   221     for (InArcIt e(gr, n); e != INVALID; ++e)

   222       sum -= flow[e];

   223     bool b = (type ==  EQ && sum == supply[n]) ||

   224              (type == GEQ && sum >= supply[n]) ||

   225              (type == LEQ && sum <= supply[n]);

   226     if (!b) return false;

   227   }

   229   return true;

   230 }

   232 // Check the feasibility of the given potentials (dual soluiton)

   233 // using the "Complementary Slackness" optimality condition

   234 template < typename GR, typename LM, typename UM,

   235            typename CM, typename SM, typename FM, typename PM >

   236 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,

   237                      const CM& cost, const SM& supply, const FM& flow,

   238                      const PM& pi, SupplyType type )

   239 {

   240   TEMPLATE_DIGRAPH_TYPEDEFS(GR);

   242   bool opt = true;

   243   for (ArcIt e(gr); opt && e != INVALID; ++e) {

   244     typename CM::Value red_cost =

   245       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];

   246     opt = red_cost == 0 ||

   247           (red_cost > 0 && flow[e] == lower[e]) ||

   248           (red_cost < 0 && flow[e] == upper[e]);

   249   }

   251   for (NodeIt n(gr); opt && n != INVALID; ++n) {

   252     typename SM::Value sum = 0;

   253     for (OutArcIt e(gr, n); e != INVALID; ++e)

   254       sum += flow[e];

   255     for (InArcIt e(gr, n); e != INVALID; ++e)

   256       sum -= flow[e];

   257     if (type != LEQ) {

   258       opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);

   259     } else {

   260       opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);

   261     }

   262   }

   264   return opt;

   265 }

   267 // Check whether the dual cost is equal to the primal cost

   268 template < typename GR, typename LM, typename UM,

   269            typename CM, typename SM, typename PM >

   270 bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,

   271                     const CM& cost, const SM& supply, const PM& pi,

   272                     typename CM::Value total )

   273 {

   274   TEMPLATE_DIGRAPH_TYPEDEFS(GR);

   276   typename CM::Value dual_cost = 0;

   277   SM red_supply(gr);

   278   for (NodeIt n(gr); n != INVALID; ++n) {

   279     red_supply[n] = supply[n];

   280   }

   281   for (ArcIt a(gr); a != INVALID; ++a) {

   282     if (lower[a] != 0) {

   283       dual_cost += lower[a] * cost[a];

   284       red_supply[gr.source(a)] -= lower[a];

   285       red_supply[gr.target(a)] += lower[a];

   286     }

   287   }

   289   for (NodeIt n(gr); n != INVALID; ++n) {

   290     dual_cost -= red_supply[n] * pi[n];

   291   }

   292   for (ArcIt a(gr); a != INVALID; ++a) {

   293     typename CM::Value red_cost =

   294       cost[a] + pi[gr.source(a)] - pi[gr.target(a)];

   295     dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);

   296   }

   298   return dual_cost == total;

   299 }

   301 // Run a minimum cost flow algorithm and check the results

   302 template < typename MCF, typename GR,

   303            typename LM, typename UM,

   304            typename CM, typename SM,

   305            typename PT >

   306 void checkMcf( const MCF& mcf, PT mcf_result,

   307                const GR& gr, const LM& lower, const UM& upper,

   308                const CM& cost, const SM& supply,

   309                PT result, bool optimal, typename CM::Value total,

   310                const std::string &test_id = "",

   311                SupplyType type = EQ )

   312 {

   313   check(mcf_result == result, "Wrong result " + test_id);

   314   if (optimal) {

   315     typename GR::template ArcMap<typename SM::Value> flow(gr);

   316     typename GR::template NodeMap<typename CM::Value> pi(gr);

   317     mcf.flowMap(flow);

   318     mcf.potentialMap(pi);

   319     check(checkFlow(gr, lower, upper, supply, flow, type),

   320           "The flow is not feasible " + test_id);

   321     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);

   322     check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),

   323           "Wrong potentials " + test_id);

   324     check(checkDualCost(gr, lower, upper, cost, supply, pi, total),

   325           "Wrong dual cost " + test_id);

   326   }

   327 }

   329 template < typename MCF, typename Param >

   330 void runMcfGeqTests( Param param,

   331                      const std::string &test_str = "",

   332                      bool full_neg_cost_support = false )

   333 {

   334   MCF mcf1(gr), mcf2(neg1_gr), mcf3(neg2_gr);

   336   // Basic tests

   337   mcf1.upperMap(u).costMap(c).supplyMap(s1);

   338   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s1,

   339            mcf1.OPTIMAL, true,     5240, test_str + "-1");

   340   mcf1.stSupply(v, w, 27);

   341   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s2,

   342            mcf1.OPTIMAL, true,     7620, test_str + "-2");

   343   mcf1.lowerMap(l2).supplyMap(s1);

   344   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s1,

   345            mcf1.OPTIMAL, true,     5970, test_str + "-3");

   346   mcf1.stSupply(v, w, 27);

   347   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s2,

   348            mcf1.OPTIMAL, true,     8010, test_str + "-4");

   349   mcf1.resetParams().supplyMap(s1);

   350   checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s1,

   351            mcf1.OPTIMAL, true,       74, test_str + "-5");

   352   mcf1.lowerMap(l2).stSupply(v, w, 27);

   353   checkMcf(mcf1, mcf1.run(param), gr, l2, cu, cc, s2,

   354            mcf1.OPTIMAL, true,       94, test_str + "-6");

   355   mcf1.reset();

   356   checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s3,

   357            mcf1.OPTIMAL, true,        0, test_str + "-7");

   358   mcf1.lowerMap(l2).upperMap(u);

   359   checkMcf(mcf1, mcf1.run(param), gr, l2, u, cc, s3,

   360            mcf1.INFEASIBLE, false,    0, test_str + "-8");

   361   mcf1.lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);

   362   checkMcf(mcf1, mcf1.run(param), gr, l3, u, c, s4,

   363            mcf1.OPTIMAL, true,     6360, test_str + "-9");

   365   // Tests for the GEQ form

   366   mcf1.resetParams().upperMap(u).costMap(c).supplyMap(s5);

   367   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s5,

   368            mcf1.OPTIMAL, true,     3530, test_str + "-10", GEQ);

   369   mcf1.lowerMap(l2);

   370   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,

   371            mcf1.OPTIMAL, true,     4540, test_str + "-11", GEQ);

   372   mcf1.supplyMap(s6);

   373   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,

   374            mcf1.INFEASIBLE, false,    0, test_str + "-12", GEQ);

   376   // Tests with negative costs

   377   mcf2.lowerMap(neg1_l1).costMap(neg1_c).supplyMap(neg1_s);

   378   checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u1, neg1_c, neg1_s,

   379            mcf2.UNBOUNDED, false,     0, test_str + "-13");

   380   mcf2.upperMap(neg1_u2);

   381   checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u2, neg1_c, neg1_s,

   382            mcf2.OPTIMAL, true,   -40000, test_str + "-14");

   383   mcf2.resetParams().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s);

   384   checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l2, neg1_u1, neg1_c, neg1_s,

   385            mcf2.UNBOUNDED, false,     0, test_str + "-15");

   387   mcf3.costMap(neg2_c).supplyMap(neg2_s);

   388   if (full_neg_cost_support) {

   389     checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,

   390              mcf3.OPTIMAL, true,   -300, test_str + "-16", GEQ);

   391   } else {

   392     checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,

   393              mcf3.UNBOUNDED, false,   0, test_str + "-17", GEQ);

   394   }

   395   mcf3.upperMap(neg2_u);

   396   checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,

   397            mcf3.OPTIMAL, true,     -300, test_str + "-18", GEQ);

   398 }

   400 template < typename MCF, typename Param >

   401 void runMcfLeqTests( Param param,

   402                      const std::string &test_str = "" )

   403 {

   404   // Tests for the LEQ form

   405   MCF mcf1(gr);

   406   mcf1.supplyType(mcf1.LEQ);

   407   mcf1.upperMap(u).costMap(c).supplyMap(s6);

   408   checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6,

   409            mcf1.OPTIMAL, true,   5080, test_str + "-19", LEQ);

   410   mcf1.lowerMap(l2);

   411   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,

   412            mcf1.OPTIMAL, true,   5930, test_str + "-20", LEQ);

   413   mcf1.supplyMap(s5);

   414   checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,

   415            mcf1.INFEASIBLE, false,  0, test_str + "-21", LEQ);

   416 }

   419 int main()

   420 {

   421   // Read the test networks

   422   std::istringstream input(test_lgf);

   423   DigraphReader<Digraph>(gr, input)

   424     .arcMap("cost", c)

   425     .arcMap("cap", u)

   426     .arcMap("low1", l1)

   427     .arcMap("low2", l2)

   428     .arcMap("low3", l3)

   429     .nodeMap("sup1", s1)

   430     .nodeMap("sup2", s2)

   431     .nodeMap("sup3", s3)

   432     .nodeMap("sup4", s4)

   433     .nodeMap("sup5", s5)

   434     .nodeMap("sup6", s6)

   435     .node("source", v)

   436     .node("target", w)

   437     .run();

   439   std::istringstream neg_inp1(test_neg1_lgf);

   440   DigraphReader<Digraph>(neg1_gr, neg_inp1)

   441     .arcMap("cost", neg1_c)

   442     .arcMap("low1", neg1_l1)

   443     .arcMap("low2", neg1_l2)

   444     .nodeMap("sup", neg1_s)

   445     .run();

   447   std::istringstream neg_inp2(test_neg2_lgf);

   448   DigraphReader<Digraph>(neg2_gr, neg_inp2)

   449     .arcMap("cost", neg2_c)

   450     .nodeMap("sup", neg2_s)

   451     .run();

   453   // Check the interface of NetworkSimplex

   454   {

   455     typedef concepts::Digraph GR;

   456     checkConcept< McfClassConcept<GR, int, int>,

   457                   NetworkSimplex<GR> >();

   458     checkConcept< McfClassConcept<GR, double, double>,

   459                   NetworkSimplex<GR, double> >();

   460     checkConcept< McfClassConcept<GR, int, double>,

   461                   NetworkSimplex<GR, int, double> >();

   462   }

   464   // Check the interface of CapacityScaling

   465   {

   466     typedef concepts::Digraph GR;

   467     checkConcept< McfClassConcept<GR, int, int>,

   468                   CapacityScaling<GR> >();

   469     checkConcept< McfClassConcept<GR, double, double>,

   470                   CapacityScaling<GR, double> >();

   471     checkConcept< McfClassConcept<GR, int, double>,

   472                   CapacityScaling<GR, int, double> >();

   473     typedef CapacityScaling<GR>::

   474       SetHeap<concepts::Heap<int, RangeMap<int> > >::Create CAS;

   475     checkConcept< McfClassConcept<GR, int, int>, CAS >();

   476   }

   478   // Check the interface of CostScaling

   479   {

   480     typedef concepts::Digraph GR;

   481     checkConcept< McfClassConcept<GR, int, int>,

   482                   CostScaling<GR> >();

   483     checkConcept< McfClassConcept<GR, double, double>,

   484                   CostScaling<GR, double> >();

   485     checkConcept< McfClassConcept<GR, int, double>,

   486                   CostScaling<GR, int, double> >();

   487     typedef CostScaling<GR>::

   488       SetLargeCost<double>::Create COS;

   489     checkConcept< McfClassConcept<GR, int, int>, COS >();

   490   }

   492   // Check the interface of CycleCanceling

   493   {

   494     typedef concepts::Digraph GR;

   495     checkConcept< McfClassConcept<GR, int, int>,

   496                   CycleCanceling<GR> >();

   497     checkConcept< McfClassConcept<GR, double, double>,

   498                   CycleCanceling<GR, double> >();

   499     checkConcept< McfClassConcept<GR, int, double>,

   500                   CycleCanceling<GR, int, double> >();

   501   }

   503   // Test NetworkSimplex

   504   {

   505     typedef NetworkSimplex<Digraph> MCF;

   506     runMcfGeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE", true);

   507     runMcfLeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE");

   508     runMcfGeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE", true);

   509     runMcfLeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE");

   510     runMcfGeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS", true);

   511     runMcfLeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS");

   512     runMcfGeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL", true);

   513     runMcfLeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL");

   514     runMcfGeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL", true);

   515     runMcfLeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL");

   516   }

   518   // Test CapacityScaling

   519   {

   520     typedef CapacityScaling<Digraph> MCF;

   521     runMcfGeqTests<MCF>(0, "SSP");

   522     runMcfGeqTests<MCF>(2, "CAS");

   523   }

   525   // Test CostScaling

   526   {

   527     typedef CostScaling<Digraph> MCF;

   528     runMcfGeqTests<MCF>(MCF::PUSH, "COS-PR");

   529     runMcfGeqTests<MCF>(MCF::AUGMENT, "COS-AR");

   530     runMcfGeqTests<MCF>(MCF::PARTIAL_AUGMENT, "COS-PAR");

   531   }

   533   // Test CycleCanceling

   534   {

   535     typedef CycleCanceling<Digraph> MCF;

   536     runMcfGeqTests<MCF>(MCF::SIMPLE_CYCLE_CANCELING, "SCC");

   537     runMcfGeqTests<MCF>(MCF::MINIMUM_MEAN_CYCLE_CANCELING, "MMCC");

   538     runMcfGeqTests<MCF>(MCF::CANCEL_AND_TIGHTEN, "CAT");

   539   }

   541   return 0;

   542 }

author	Peter Kovacs <kpeter@inf.elte.hu>
	Sat, 16 Mar 2013 14:11:32 +0100
changeset 1220	eb2f9d453070
parent 898	75c97c3786d6
child 1270	dceba191c00d
child 1317	b40c2bbb8da5
permissions	-rw-r--r--