COIN-OR::LEMON - Graph Library

source: lemon/test/min_cost_flow_test.cc @ 689:111698359429

Last change on this file since 689:111698359429 was 689:111698359429, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Less map copying in NetworkSimplex? (#234)

  • The graph is copied in the constructor instead of the init() function. It must not be modified after the class is constructed.
  • The maps are copied once (instead of twice).
  • Remove FlowMap?, PotentialMap? typedefs and flowMap(), pontentialMap() setter functions.
  • flowMap() and potentialMap() query functions copy the values into the given map (reference) instead of returning a const reference to a previously constructed map.
File size: 12.1 KB
Line 
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
33using namespace lemon;
34
35char 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
80enum SupplyType {
81  EQ,
82  GEQ,
83  LEQ
84};
85
86// Check the interface of an MCF algorithm
87template <typename GR, typename Value, typename Cost>
88class McfClassConcept
89{
90public:
91
92  template <typename MCF>
93  struct Constraints {
94    void constraints() {
95      checkConcept<concepts::Digraph, GR>();
96
97      MCF mcf(g);
98      const MCF& const_mcf = mcf;
99
100      b = mcf.reset()
101             .lowerMap(lower)
102             .upperMap(upper)
103             .costMap(cost)
104             .supplyMap(sup)
105             .stSupply(n, n, k)
106             .run();
107
108      c = const_mcf.totalCost();
109      x = const_mcf.template totalCost<double>();
110      v = const_mcf.flow(a);
111      c = const_mcf.potential(n);
112      const_mcf.flowMap(fm);
113      const_mcf.potentialMap(pm);
114    }
115
116    typedef typename GR::Node Node;
117    typedef typename GR::Arc Arc;
118    typedef concepts::ReadMap<Node, Value> NM;
119    typedef concepts::ReadMap<Arc, Value> VAM;
120    typedef concepts::ReadMap<Arc, Cost> CAM;
121    typedef concepts::WriteMap<Arc, Value> FlowMap;
122    typedef concepts::WriteMap<Node, Cost> PotMap;
123
124    const GR &g;
125    const VAM &lower;
126    const VAM &upper;
127    const CAM &cost;
128    const NM &sup;
129    const Node &n;
130    const Arc &a;
131    const Value &k;
132    FlowMap fm;
133    PotMap pm;
134    bool b;
135    double x;
136    typename MCF::Value v;
137    typename MCF::Cost c;
138  };
139
140};
141
142
143// Check the feasibility of the given flow (primal soluiton)
144template < typename GR, typename LM, typename UM,
145           typename SM, typename FM >
146bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
147                const SM& supply, const FM& flow,
148                SupplyType type = EQ )
149{
150  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
151
152  for (ArcIt e(gr); e != INVALID; ++e) {
153    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
154  }
155
156  for (NodeIt n(gr); n != INVALID; ++n) {
157    typename SM::Value sum = 0;
158    for (OutArcIt e(gr, n); e != INVALID; ++e)
159      sum += flow[e];
160    for (InArcIt e(gr, n); e != INVALID; ++e)
161      sum -= flow[e];
162    bool b = (type ==  EQ && sum == supply[n]) ||
163             (type == GEQ && sum >= supply[n]) ||
164             (type == LEQ && sum <= supply[n]);
165    if (!b) return false;
166  }
167
168  return true;
169}
170
171// Check the feasibility of the given potentials (dual soluiton)
172// using the "Complementary Slackness" optimality condition
173template < typename GR, typename LM, typename UM,
174           typename CM, typename SM, typename FM, typename PM >
175bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
176                     const CM& cost, const SM& supply, const FM& flow,
177                     const PM& pi )
178{
179  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
180
181  bool opt = true;
182  for (ArcIt e(gr); opt && e != INVALID; ++e) {
183    typename CM::Value red_cost =
184      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
185    opt = red_cost == 0 ||
186          (red_cost > 0 && flow[e] == lower[e]) ||
187          (red_cost < 0 && flow[e] == upper[e]);
188  }
189 
190  for (NodeIt n(gr); opt && n != INVALID; ++n) {
191    typename SM::Value sum = 0;
192    for (OutArcIt e(gr, n); e != INVALID; ++e)
193      sum += flow[e];
194    for (InArcIt e(gr, n); e != INVALID; ++e)
195      sum -= flow[e];
196    opt = (sum == supply[n]) || (pi[n] == 0);
197  }
198 
199  return opt;
200}
201
202// Run a minimum cost flow algorithm and check the results
203template < typename MCF, typename GR,
204           typename LM, typename UM,
205           typename CM, typename SM,
206           typename PT >
207void checkMcf( const MCF& mcf, PT mcf_result,
208               const GR& gr, const LM& lower, const UM& upper,
209               const CM& cost, const SM& supply,
210               PT result, bool optimal, typename CM::Value total,
211               const std::string &test_id = "",
212               SupplyType type = EQ )
213{
214  check(mcf_result == result, "Wrong result " + test_id);
215  if (optimal) {
216    typename GR::template ArcMap<typename SM::Value> flow(gr);
217    typename GR::template NodeMap<typename CM::Value> pi(gr);
218    mcf.flowMap(flow);
219    mcf.potentialMap(pi);
220    check(checkFlow(gr, lower, upper, supply, flow, type),
221          "The flow is not feasible " + test_id);
222    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
223    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
224          "Wrong potentials " + test_id);
225  }
226}
227
228int main()
229{
230  // Check the interfaces
231  {
232    typedef concepts::Digraph GR;
233    checkConcept< McfClassConcept<GR, int, int>,
234                  NetworkSimplex<GR> >();
235    checkConcept< McfClassConcept<GR, double, double>,
236                  NetworkSimplex<GR, double> >();
237    checkConcept< McfClassConcept<GR, int, double>,
238                  NetworkSimplex<GR, int, double> >();
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), l3(gr), u(gr);
248  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(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    .arcMap("low3", l3)
259    .nodeMap("sup1", s1)
260    .nodeMap("sup2", s2)
261    .nodeMap("sup3", s3)
262    .nodeMap("sup4", s4)
263    .nodeMap("sup5", s5)
264    .nodeMap("sup6", s6)
265    .node("source", v)
266    .node("target", w)
267    .run();
268 
269  // Build a test digraph for testing negative costs
270  Digraph ngr;
271  Node n1 = ngr.addNode();
272  Node n2 = ngr.addNode();
273  Node n3 = ngr.addNode();
274  Node n4 = ngr.addNode();
275  Node n5 = ngr.addNode();
276  Node n6 = ngr.addNode();
277  Node n7 = ngr.addNode();
278 
279  Arc a1 = ngr.addArc(n1, n2);
280  Arc a2 = ngr.addArc(n1, n3);
281  Arc a3 = ngr.addArc(n2, n4);
282  Arc a4 = ngr.addArc(n3, n4);
283  Arc a5 = ngr.addArc(n3, n2);
284  Arc a6 = ngr.addArc(n5, n3);
285  Arc a7 = ngr.addArc(n5, n6);
286  Arc a8 = ngr.addArc(n6, n7);
287  Arc a9 = ngr.addArc(n7, n5);
288 
289  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
290  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
291  Digraph::NodeMap<int> ns(ngr, 0);
292 
293  nl2[a7] =  1000;
294  nl2[a8] = -1000;
295 
296  ns[n1] =  100;
297  ns[n4] = -100;
298 
299  nc[a1] =  100;
300  nc[a2] =   30;
301  nc[a3] =   20;
302  nc[a4] =   80;
303  nc[a5] =   50;
304  nc[a6] =   10;
305  nc[a7] =   80;
306  nc[a8] =   30;
307  nc[a9] = -120;
308
309  // A. Test NetworkSimplex with the default pivot rule
310  {
311    NetworkSimplex<Digraph> mcf(gr);
312
313    // Check the equality form
314    mcf.upperMap(u).costMap(c);
315    checkMcf(mcf, mcf.supplyMap(s1).run(),
316             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
317    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
318             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
319    mcf.lowerMap(l2);
320    checkMcf(mcf, mcf.supplyMap(s1).run(),
321             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
322    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
323             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
324    mcf.reset();
325    checkMcf(mcf, mcf.supplyMap(s1).run(),
326             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
327    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
328             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
329    mcf.reset();
330    checkMcf(mcf, mcf.run(),
331             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
332    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
333             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
334    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
335    checkMcf(mcf, mcf.run(),
336             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
337
338    // Check the GEQ form
339    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
340    checkMcf(mcf, mcf.run(),
341             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
342    mcf.supplyType(mcf.GEQ);
343    checkMcf(mcf, mcf.lowerMap(l2).run(),
344             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
345    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
346    checkMcf(mcf, mcf.run(),
347             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
348
349    // Check the LEQ form
350    mcf.reset().supplyType(mcf.LEQ);
351    mcf.upperMap(u).costMap(c).supplyMap(s6);
352    checkMcf(mcf, mcf.run(),
353             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
354    checkMcf(mcf, mcf.lowerMap(l2).run(),
355             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
356    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
357    checkMcf(mcf, mcf.run(),
358             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
359
360    // Check negative costs
361    NetworkSimplex<Digraph> nmcf(ngr);
362    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
363    checkMcf(nmcf, nmcf.run(),
364      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
365    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
366      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
367    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
368    checkMcf(nmcf, nmcf.run(),
369      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
370  }
371
372  // B. Test NetworkSimplex with each pivot rule
373  {
374    NetworkSimplex<Digraph> mcf(gr);
375    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
376
377    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
378             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
379    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
380             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
381    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
382             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
383    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
384             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
385    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
386             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
387  }
388
389  return 0;
390}
Note: See TracBrowser for help on using the repository browser.