COIN-OR::LEMON - Graph Library

source: lemon/test/min_cost_flow_test.cc @ 658:85cb3aa71cce

Last change on this file since 658:85cb3aa71cce was 656:e6927fe719e6, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Support >= and <= constraints in NetworkSimplex? (#219, #234)

By default the same inequality constraints are supported as by
Circulation (the GEQ form), but the LEQ form can also be selected
using the problemType() function.

The documentation of the min. cost flow module is reworked and
extended with important notes and explanations about the different
variants of the problem and about the dual solution and optimality
conditions.

File size: 10.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
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
32using namespace lemon;
33
34char 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
79enum ProblemType {
80  EQ,
81  GEQ,
82  LEQ
83};
84
85// Check the interface of an MCF algorithm
86template <typename GR, typename Flow, typename Cost>
87class McfClassConcept
88{
89public:
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)
150template < typename GR, typename LM, typename UM,
151           typename SM, typename FM >
152bool 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
179template < typename GR, typename LM, typename UM,
180           typename CM, typename SM, typename FM, typename PM >
181bool 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
209template < typename MCF, typename GR,
210           typename LM, typename UM,
211           typename CM, typename SM >
212void 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
230int 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}
Note: See TracBrowser for help on using the repository browser.