COIN-OR::LEMON - Graph Library

source: lemon-1.2/test/min_cost_flow_test.cc @ 641:756a5ec551c8

Last change on this file since 641:756a5ec551c8 was 640:6c408d864fa1, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Support negative costs and bounds in NetworkSimplex? (#270)

  • The interface is reworked to support negative costs and bounds.
    • ProblemType? and problemType() are renamed to SupplyType? and supplyType(), see also #234.
    • ProblemType? type is introduced similarly to the LP interface.
    • 'bool run()' is replaced by 'ProblemType? run()' to handle unbounded problem instances, as well.
    • Add INF public member constant similarly to the LP interface.
  • Remove capacityMap() and boundMaps(), see also #266.
  • Update the problem definition in the MCF module.
  • Remove the usage of Circulation (and adaptors) for checking feasibility. Check feasibility by examining the artifical arcs instead (after solving the problem).
  • Additional check for unbounded negative cycles found during the algorithm (it is possible now, since negative costs are allowed).
  • Fix in the constructor (the value types needn't be integer any more), see also #254.
  • Improve and extend the doc.
  • Rework the test file and add test cases for negative costs and bounds.
File size: 12.0 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 Flow, 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
99      b = mcf.reset()
100             .lowerMap(lower)
101             .upperMap(upper)
102             .costMap(cost)
103             .supplyMap(sup)
104             .stSupply(n, n, k)
105             .flowMap(flow)
106             .potentialMap(pot)
107             .run();
108     
109      const MCF& const_mcf = mcf;
110
111      const typename MCF::FlowMap &fm = const_mcf.flowMap();
112      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
113
114      c = const_mcf.totalCost();
115      double x = const_mcf.template totalCost<double>();
116      v = const_mcf.flow(a);
117      c = const_mcf.potential(n);
118     
119      v = const_mcf.INF;
120
121      ignore_unused_variable_warning(fm);
122      ignore_unused_variable_warning(pm);
123      ignore_unused_variable_warning(x);
124    }
125
126    typedef typename GR::Node Node;
127    typedef typename GR::Arc Arc;
128    typedef concepts::ReadMap<Node, Flow> NM;
129    typedef concepts::ReadMap<Arc, Flow> FAM;
130    typedef concepts::ReadMap<Arc, Cost> CAM;
131
132    const GR &g;
133    const FAM &lower;
134    const FAM &upper;
135    const CAM &cost;
136    const NM &sup;
137    const Node &n;
138    const Arc &a;
139    const Flow &k;
140    Flow v;
141    Cost c;
142    bool b;
143
144    typename MCF::FlowMap &flow;
145    typename MCF::PotentialMap &pot;
146  };
147
148};
149
150
151// Check the feasibility of the given flow (primal soluiton)
152template < typename GR, typename LM, typename UM,
153           typename SM, typename FM >
154bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
155                const SM& supply, const FM& flow,
156                SupplyType type = EQ )
157{
158  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
159
160  for (ArcIt e(gr); e != INVALID; ++e) {
161    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
162  }
163
164  for (NodeIt n(gr); n != INVALID; ++n) {
165    typename SM::Value sum = 0;
166    for (OutArcIt e(gr, n); e != INVALID; ++e)
167      sum += flow[e];
168    for (InArcIt e(gr, n); e != INVALID; ++e)
169      sum -= flow[e];
170    bool b = (type ==  EQ && sum == supply[n]) ||
171             (type == GEQ && sum >= supply[n]) ||
172             (type == LEQ && sum <= supply[n]);
173    if (!b) return false;
174  }
175
176  return true;
177}
178
179// Check the feasibility of the given potentials (dual soluiton)
180// using the "Complementary Slackness" optimality condition
181template < typename GR, typename LM, typename UM,
182           typename CM, typename SM, typename FM, typename PM >
183bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
184                     const CM& cost, const SM& supply, const FM& flow,
185                     const PM& pi )
186{
187  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
188
189  bool opt = true;
190  for (ArcIt e(gr); opt && e != INVALID; ++e) {
191    typename CM::Value red_cost =
192      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
193    opt = red_cost == 0 ||
194          (red_cost > 0 && flow[e] == lower[e]) ||
195          (red_cost < 0 && flow[e] == upper[e]);
196  }
197 
198  for (NodeIt n(gr); opt && n != INVALID; ++n) {
199    typename SM::Value sum = 0;
200    for (OutArcIt e(gr, n); e != INVALID; ++e)
201      sum += flow[e];
202    for (InArcIt e(gr, n); e != INVALID; ++e)
203      sum -= flow[e];
204    opt = (sum == supply[n]) || (pi[n] == 0);
205  }
206 
207  return opt;
208}
209
210// Run a minimum cost flow algorithm and check the results
211template < typename MCF, typename GR,
212           typename LM, typename UM,
213           typename CM, typename SM,
214           typename PT >
215void checkMcf( const MCF& mcf, PT mcf_result,
216               const GR& gr, const LM& lower, const UM& upper,
217               const CM& cost, const SM& supply,
218               PT result, bool optimal, typename CM::Value total,
219               const std::string &test_id = "",
220               SupplyType type = EQ )
221{
222  check(mcf_result == result, "Wrong result " + test_id);
223  if (optimal) {
224    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
225          "The flow is not feasible " + test_id);
226    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
227    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
228                         mcf.potentialMap()),
229          "Wrong potentials " + test_id);
230  }
231}
232
233int main()
234{
235  // Check the interfaces
236  {
237    typedef int Flow;
238    typedef int Cost;
239    typedef concepts::Digraph GR;
240    checkConcept< McfClassConcept<GR, Flow, Cost>,
241                  NetworkSimplex<GR, Flow, Cost> >();
242  }
243
244  // Run various MCF tests
245  typedef ListDigraph Digraph;
246  DIGRAPH_TYPEDEFS(ListDigraph);
247
248  // Read the test digraph
249  Digraph gr;
250  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
251  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
252  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
253  Node v, w;
254
255  std::istringstream input(test_lgf);
256  DigraphReader<Digraph>(gr, input)
257    .arcMap("cost", c)
258    .arcMap("cap", u)
259    .arcMap("low1", l1)
260    .arcMap("low2", l2)
261    .arcMap("low3", l3)
262    .nodeMap("sup1", s1)
263    .nodeMap("sup2", s2)
264    .nodeMap("sup3", s3)
265    .nodeMap("sup4", s4)
266    .nodeMap("sup5", s5)
267    .nodeMap("sup6", s6)
268    .node("source", v)
269    .node("target", w)
270    .run();
271 
272  // Build a test digraph for testing negative costs
273  Digraph ngr;
274  Node n1 = ngr.addNode();
275  Node n2 = ngr.addNode();
276  Node n3 = ngr.addNode();
277  Node n4 = ngr.addNode();
278  Node n5 = ngr.addNode();
279  Node n6 = ngr.addNode();
280  Node n7 = ngr.addNode();
281 
282  Arc a1 = ngr.addArc(n1, n2);
283  Arc a2 = ngr.addArc(n1, n3);
284  Arc a3 = ngr.addArc(n2, n4);
285  Arc a4 = ngr.addArc(n3, n4);
286  Arc a5 = ngr.addArc(n3, n2);
287  Arc a6 = ngr.addArc(n5, n3);
288  Arc a7 = ngr.addArc(n5, n6);
289  Arc a8 = ngr.addArc(n6, n7);
290  Arc a9 = ngr.addArc(n7, n5);
291 
292  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
293  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
294  Digraph::NodeMap<int> ns(ngr, 0);
295 
296  nl2[a7] =  1000;
297  nl2[a8] = -1000;
298 
299  ns[n1] =  100;
300  ns[n4] = -100;
301 
302  nc[a1] =  100;
303  nc[a2] =   30;
304  nc[a3] =   20;
305  nc[a4] =   80;
306  nc[a5] =   50;
307  nc[a6] =   10;
308  nc[a7] =   80;
309  nc[a8] =   30;
310  nc[a9] = -120;
311
312  // A. Test NetworkSimplex with the default pivot rule
313  {
314    NetworkSimplex<Digraph> mcf(gr);
315
316    // Check the equality form
317    mcf.upperMap(u).costMap(c);
318    checkMcf(mcf, mcf.supplyMap(s1).run(),
319             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
320    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
321             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
322    mcf.lowerMap(l2);
323    checkMcf(mcf, mcf.supplyMap(s1).run(),
324             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
325    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
326             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
327    mcf.reset();
328    checkMcf(mcf, mcf.supplyMap(s1).run(),
329             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
330    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
331             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
332    mcf.reset();
333    checkMcf(mcf, mcf.run(),
334             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
335    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
336             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
337    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
338    checkMcf(mcf, mcf.run(),
339             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
340
341    // Check the GEQ form
342    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
343    checkMcf(mcf, mcf.run(),
344             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
345    mcf.supplyType(mcf.GEQ);
346    checkMcf(mcf, mcf.lowerMap(l2).run(),
347             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
348    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
349    checkMcf(mcf, mcf.run(),
350             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
351
352    // Check the LEQ form
353    mcf.reset().supplyType(mcf.LEQ);
354    mcf.upperMap(u).costMap(c).supplyMap(s6);
355    checkMcf(mcf, mcf.run(),
356             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
357    checkMcf(mcf, mcf.lowerMap(l2).run(),
358             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
359    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
360    checkMcf(mcf, mcf.run(),
361             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
362
363    // Check negative costs
364    NetworkSimplex<Digraph> nmcf(ngr);
365    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
366    checkMcf(nmcf, nmcf.run(),
367      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
368    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
369      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
370    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
371    checkMcf(nmcf, nmcf.run(),
372      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
373  }
374
375  // B. Test NetworkSimplex with each pivot rule
376  {
377    NetworkSimplex<Digraph> mcf(gr);
378    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
379
380    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
381             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
382    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
383             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
384    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
385             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
386    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
387             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
388    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
389             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
390  }
391
392  return 0;
393}
Note: See TracBrowser for help on using the repository browser.