COIN-OR::LEMON - Graph Library

source: lemon-1.2/test/min_cost_flow_test.cc @ 781:6f10c6ec5a21

Last change on this file since 781:6f10c6ec5a21 was 669:4faca85d40e6, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Avoid Intel C++ Compiler warnings

File size: 14.0 KB
RevLine 
[601]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>
[640]21#include <limits>
[601]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
33using namespace lemon;
34
35char test_lgf[] =
36  "@nodes\n"
[640]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"               
[601]51  "@arcs\n"
[640]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"
[601]74  "\n"
75  "@attributes\n"
76  "source 1\n"
77  "target 12\n";
78
79
[640]80enum SupplyType {
[609]81  EQ,
82  GEQ,
83  LEQ
84};
85
[601]86// Check the interface of an MCF algorithm
[642]87template <typename GR, typename Value, typename Cost>
[601]88class McfClassConcept
89{
90public:
91
92  template <typename MCF>
93  struct Constraints {
94    void constraints() {
95      checkConcept<concepts::Digraph, GR>();
[669]96     
97      const Constraints& me = *this;
[601]98
[669]99      MCF mcf(me.g);
[642]100      const MCF& const_mcf = mcf;
[601]101
[606]102      b = mcf.reset()
[669]103             .lowerMap(me.lower)
104             .upperMap(me.upper)
105             .costMap(me.cost)
106             .supplyMap(me.sup)
107             .stSupply(me.n, me.n, me.k)
[605]108             .run();
109
[640]110      c = const_mcf.totalCost();
[642]111      x = const_mcf.template totalCost<double>();
[669]112      v = const_mcf.flow(me.a);
113      c = const_mcf.potential(me.n);
[642]114      const_mcf.flowMap(fm);
115      const_mcf.potentialMap(pm);
[601]116    }
117
118    typedef typename GR::Node Node;
119    typedef typename GR::Arc Arc;
[642]120    typedef concepts::ReadMap<Node, Value> NM;
121    typedef concepts::ReadMap<Arc, Value> VAM;
[607]122    typedef concepts::ReadMap<Arc, Cost> CAM;
[642]123    typedef concepts::WriteMap<Arc, Value> FlowMap;
124    typedef concepts::WriteMap<Node, Cost> PotMap;
[669]125 
126    GR g;
127    VAM lower;
128    VAM upper;
129    CAM cost;
130    NM sup;
131    Node n;
132    Arc a;
133    Value k;
[601]134
[642]135    FlowMap fm;
136    PotMap pm;
[605]137    bool b;
[642]138    double x;
139    typename MCF::Value v;
140    typename MCF::Cost c;
[601]141  };
142
143};
144
145
146// Check the feasibility of the given flow (primal soluiton)
147template < typename GR, typename LM, typename UM,
148           typename SM, typename FM >
149bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
[609]150                const SM& supply, const FM& flow,
[640]151                SupplyType type = EQ )
[601]152{
153  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
154
155  for (ArcIt e(gr); e != INVALID; ++e) {
156    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
157  }
158
159  for (NodeIt n(gr); n != INVALID; ++n) {
160    typename SM::Value sum = 0;
161    for (OutArcIt e(gr, n); e != INVALID; ++e)
162      sum += flow[e];
163    for (InArcIt e(gr, n); e != INVALID; ++e)
164      sum -= flow[e];
[609]165    bool b = (type ==  EQ && sum == supply[n]) ||
166             (type == GEQ && sum >= supply[n]) ||
167             (type == LEQ && sum <= supply[n]);
168    if (!b) return false;
[601]169  }
170
171  return true;
172}
173
174// Check the feasibility of the given potentials (dual soluiton)
[605]175// using the "Complementary Slackness" optimality condition
[601]176template < typename GR, typename LM, typename UM,
[609]177           typename CM, typename SM, typename FM, typename PM >
[601]178bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
[609]179                     const CM& cost, const SM& supply, const FM& flow,
[664]180                     const PM& pi, SupplyType type )
[601]181{
182  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
183
184  bool opt = true;
185  for (ArcIt e(gr); opt && e != INVALID; ++e) {
186    typename CM::Value red_cost =
187      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
188    opt = red_cost == 0 ||
189          (red_cost > 0 && flow[e] == lower[e]) ||
190          (red_cost < 0 && flow[e] == upper[e]);
191  }
[609]192 
193  for (NodeIt n(gr); opt && n != INVALID; ++n) {
194    typename SM::Value sum = 0;
195    for (OutArcIt e(gr, n); e != INVALID; ++e)
196      sum += flow[e];
197    for (InArcIt e(gr, n); e != INVALID; ++e)
198      sum -= flow[e];
[664]199    if (type != LEQ) {
200      opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
201    } else {
202      opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
203    }
[609]204  }
205 
[601]206  return opt;
207}
208
[664]209// Check whether the dual cost is equal to the primal cost
210template < typename GR, typename LM, typename UM,
211           typename CM, typename SM, typename PM >
212bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
213                    const CM& cost, const SM& supply, const PM& pi,
214                    typename CM::Value total )
215{
216  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
217
218  typename CM::Value dual_cost = 0;
219  SM red_supply(gr);
220  for (NodeIt n(gr); n != INVALID; ++n) {
221    red_supply[n] = supply[n];
222  }
223  for (ArcIt a(gr); a != INVALID; ++a) {
224    if (lower[a] != 0) {
225      dual_cost += lower[a] * cost[a];
226      red_supply[gr.source(a)] -= lower[a];
227      red_supply[gr.target(a)] += lower[a];
228    }
229  }
230 
231  for (NodeIt n(gr); n != INVALID; ++n) {
232    dual_cost -= red_supply[n] * pi[n];
233  }
234  for (ArcIt a(gr); a != INVALID; ++a) {
235    typename CM::Value red_cost =
236      cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
237    dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
238  }
239 
240  return dual_cost == total;
241}
242
[601]243// Run a minimum cost flow algorithm and check the results
244template < typename MCF, typename GR,
245           typename LM, typename UM,
[640]246           typename CM, typename SM,
247           typename PT >
248void checkMcf( const MCF& mcf, PT mcf_result,
[601]249               const GR& gr, const LM& lower, const UM& upper,
250               const CM& cost, const SM& supply,
[640]251               PT result, bool optimal, typename CM::Value total,
[609]252               const std::string &test_id = "",
[640]253               SupplyType type = EQ )
[601]254{
255  check(mcf_result == result, "Wrong result " + test_id);
[640]256  if (optimal) {
[642]257    typename GR::template ArcMap<typename SM::Value> flow(gr);
258    typename GR::template NodeMap<typename CM::Value> pi(gr);
259    mcf.flowMap(flow);
260    mcf.potentialMap(pi);
261    check(checkFlow(gr, lower, upper, supply, flow, type),
[601]262          "The flow is not feasible " + test_id);
263    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
[664]264    check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
[601]265          "Wrong potentials " + test_id);
[664]266    check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
267          "Wrong dual cost " + test_id);
[601]268  }
269}
270
271int main()
272{
273  // Check the interfaces
274  {
[615]275    typedef concepts::Digraph GR;
[642]276    checkConcept< McfClassConcept<GR, int, int>,
277                  NetworkSimplex<GR> >();
278    checkConcept< McfClassConcept<GR, double, double>,
279                  NetworkSimplex<GR, double> >();
280    checkConcept< McfClassConcept<GR, int, double>,
281                  NetworkSimplex<GR, int, double> >();
[601]282  }
283
284  // Run various MCF tests
285  typedef ListDigraph Digraph;
286  DIGRAPH_TYPEDEFS(ListDigraph);
287
288  // Read the test digraph
289  Digraph gr;
[640]290  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
291  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
[605]292  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
[601]293  Node v, w;
294
295  std::istringstream input(test_lgf);
296  DigraphReader<Digraph>(gr, input)
297    .arcMap("cost", c)
298    .arcMap("cap", u)
299    .arcMap("low1", l1)
300    .arcMap("low2", l2)
[640]301    .arcMap("low3", l3)
[601]302    .nodeMap("sup1", s1)
303    .nodeMap("sup2", s2)
304    .nodeMap("sup3", s3)
[609]305    .nodeMap("sup4", s4)
306    .nodeMap("sup5", s5)
[640]307    .nodeMap("sup6", s6)
[601]308    .node("source", v)
309    .node("target", w)
310    .run();
[640]311 
[664]312  // Build test digraphs with negative costs
313  Digraph neg_gr;
314  Node n1 = neg_gr.addNode();
315  Node n2 = neg_gr.addNode();
316  Node n3 = neg_gr.addNode();
317  Node n4 = neg_gr.addNode();
318  Node n5 = neg_gr.addNode();
319  Node n6 = neg_gr.addNode();
320  Node n7 = neg_gr.addNode();
[640]321 
[664]322  Arc a1 = neg_gr.addArc(n1, n2);
323  Arc a2 = neg_gr.addArc(n1, n3);
324  Arc a3 = neg_gr.addArc(n2, n4);
325  Arc a4 = neg_gr.addArc(n3, n4);
326  Arc a5 = neg_gr.addArc(n3, n2);
327  Arc a6 = neg_gr.addArc(n5, n3);
328  Arc a7 = neg_gr.addArc(n5, n6);
329  Arc a8 = neg_gr.addArc(n6, n7);
330  Arc a9 = neg_gr.addArc(n7, n5);
[640]331 
[664]332  Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
333  ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
334  Digraph::NodeMap<int> neg_s(neg_gr, 0);
[640]335 
[664]336  neg_l2[a7] =  1000;
337  neg_l2[a8] = -1000;
[640]338 
[664]339  neg_s[n1] =  100;
340  neg_s[n4] = -100;
[640]341 
[664]342  neg_c[a1] =  100;
343  neg_c[a2] =   30;
344  neg_c[a3] =   20;
345  neg_c[a4] =   80;
346  neg_c[a5] =   50;
347  neg_c[a6] =   10;
348  neg_c[a7] =   80;
349  neg_c[a8] =   30;
350  neg_c[a9] = -120;
351
352  Digraph negs_gr;
353  Digraph::NodeMap<int> negs_s(negs_gr);
354  Digraph::ArcMap<int> negs_c(negs_gr);
355  ConstMap<Arc, int> negs_l(0), negs_u(1000);
356  n1 = negs_gr.addNode();
357  n2 = negs_gr.addNode();
358  negs_s[n1] = 100;
359  negs_s[n2] = -300;
360  negs_c[negs_gr.addArc(n1, n2)] = -1;
361
[601]362
[605]363  // A. Test NetworkSimplex with the default pivot rule
[601]364  {
[606]365    NetworkSimplex<Digraph> mcf(gr);
[601]366
[609]367    // Check the equality form
[606]368    mcf.upperMap(u).costMap(c);
369    checkMcf(mcf, mcf.supplyMap(s1).run(),
[640]370             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
[606]371    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
[640]372             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
[606]373    mcf.lowerMap(l2);
374    checkMcf(mcf, mcf.supplyMap(s1).run(),
[640]375             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
[606]376    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
[640]377             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
[606]378    mcf.reset();
379    checkMcf(mcf, mcf.supplyMap(s1).run(),
[640]380             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
[606]381    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
[640]382             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
[606]383    mcf.reset();
384    checkMcf(mcf, mcf.run(),
[640]385             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
386    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
387             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
388    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
389    checkMcf(mcf, mcf.run(),
390             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
[609]391
392    // Check the GEQ form
[640]393    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
[609]394    checkMcf(mcf, mcf.run(),
[640]395             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
396    mcf.supplyType(mcf.GEQ);
[609]397    checkMcf(mcf, mcf.lowerMap(l2).run(),
[640]398             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
[664]399    mcf.supplyMap(s6);
[609]400    checkMcf(mcf, mcf.run(),
[640]401             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
[609]402
403    // Check the LEQ form
[640]404    mcf.reset().supplyType(mcf.LEQ);
405    mcf.upperMap(u).costMap(c).supplyMap(s6);
[609]406    checkMcf(mcf, mcf.run(),
[640]407             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
[609]408    checkMcf(mcf, mcf.lowerMap(l2).run(),
[640]409             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
[664]410    mcf.supplyMap(s5);
[609]411    checkMcf(mcf, mcf.run(),
[640]412             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
413
414    // Check negative costs
[664]415    NetworkSimplex<Digraph> neg_mcf(neg_gr);
416    neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
417    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
418      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
419    neg_mcf.upperMap(neg_u2);
420    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
421      neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
422    neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
423    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
424      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");
425     
426    NetworkSimplex<Digraph> negs_mcf(negs_gr);
427    negs_mcf.costMap(negs_c).supplyMap(negs_s);
428    checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
429      negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);
[601]430  }
431
[605]432  // B. Test NetworkSimplex with each pivot rule
[601]433  {
[606]434    NetworkSimplex<Digraph> mcf(gr);
[640]435    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
[601]436
[606]437    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
[640]438             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
[606]439    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
[640]440             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
[606]441    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
[640]442             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
[606]443    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
[640]444             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
[606]445    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
[640]446             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
[601]447  }
448
449  return 0;
450}
Note: See TracBrowser for help on using the repository browser.