COIN-OR::LEMON - Graph Library

source: lemon-1.2/test/min_cost_flow_test.cc @ 668:7e13120d90a2

Last change on this file since 668:7e13120d90a2 was 664:cc61d09f053b, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Extend min cost flow test file + check dual costs (#291)

File size: 14.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 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, SupplyType type )
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    if (type != LEQ) {
197      opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
198    } else {
199      opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
200    }
201  }
202 
203  return opt;
204}
205
206// Check whether the dual cost is equal to the primal cost
207template < typename GR, typename LM, typename UM,
208           typename CM, typename SM, typename PM >
209bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
210                    const CM& cost, const SM& supply, const PM& pi,
211                    typename CM::Value total )
212{
213  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
214
215  typename CM::Value dual_cost = 0;
216  SM red_supply(gr);
217  for (NodeIt n(gr); n != INVALID; ++n) {
218    red_supply[n] = supply[n];
219  }
220  for (ArcIt a(gr); a != INVALID; ++a) {
221    if (lower[a] != 0) {
222      dual_cost += lower[a] * cost[a];
223      red_supply[gr.source(a)] -= lower[a];
224      red_supply[gr.target(a)] += lower[a];
225    }
226  }
227 
228  for (NodeIt n(gr); n != INVALID; ++n) {
229    dual_cost -= red_supply[n] * pi[n];
230  }
231  for (ArcIt a(gr); a != INVALID; ++a) {
232    typename CM::Value red_cost =
233      cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
234    dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
235  }
236 
237  return dual_cost == total;
238}
239
240// Run a minimum cost flow algorithm and check the results
241template < typename MCF, typename GR,
242           typename LM, typename UM,
243           typename CM, typename SM,
244           typename PT >
245void checkMcf( const MCF& mcf, PT mcf_result,
246               const GR& gr, const LM& lower, const UM& upper,
247               const CM& cost, const SM& supply,
248               PT result, bool optimal, typename CM::Value total,
249               const std::string &test_id = "",
250               SupplyType type = EQ )
251{
252  check(mcf_result == result, "Wrong result " + test_id);
253  if (optimal) {
254    typename GR::template ArcMap<typename SM::Value> flow(gr);
255    typename GR::template NodeMap<typename CM::Value> pi(gr);
256    mcf.flowMap(flow);
257    mcf.potentialMap(pi);
258    check(checkFlow(gr, lower, upper, supply, flow, type),
259          "The flow is not feasible " + test_id);
260    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
261    check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
262          "Wrong potentials " + test_id);
263    check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
264          "Wrong dual cost " + test_id);
265  }
266}
267
268int main()
269{
270  // Check the interfaces
271  {
272    typedef concepts::Digraph GR;
273    checkConcept< McfClassConcept<GR, int, int>,
274                  NetworkSimplex<GR> >();
275    checkConcept< McfClassConcept<GR, double, double>,
276                  NetworkSimplex<GR, double> >();
277    checkConcept< McfClassConcept<GR, int, double>,
278                  NetworkSimplex<GR, int, double> >();
279  }
280
281  // Run various MCF tests
282  typedef ListDigraph Digraph;
283  DIGRAPH_TYPEDEFS(ListDigraph);
284
285  // Read the test digraph
286  Digraph gr;
287  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
288  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
289  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
290  Node v, w;
291
292  std::istringstream input(test_lgf);
293  DigraphReader<Digraph>(gr, input)
294    .arcMap("cost", c)
295    .arcMap("cap", u)
296    .arcMap("low1", l1)
297    .arcMap("low2", l2)
298    .arcMap("low3", l3)
299    .nodeMap("sup1", s1)
300    .nodeMap("sup2", s2)
301    .nodeMap("sup3", s3)
302    .nodeMap("sup4", s4)
303    .nodeMap("sup5", s5)
304    .nodeMap("sup6", s6)
305    .node("source", v)
306    .node("target", w)
307    .run();
308 
309  // Build test digraphs with negative costs
310  Digraph neg_gr;
311  Node n1 = neg_gr.addNode();
312  Node n2 = neg_gr.addNode();
313  Node n3 = neg_gr.addNode();
314  Node n4 = neg_gr.addNode();
315  Node n5 = neg_gr.addNode();
316  Node n6 = neg_gr.addNode();
317  Node n7 = neg_gr.addNode();
318 
319  Arc a1 = neg_gr.addArc(n1, n2);
320  Arc a2 = neg_gr.addArc(n1, n3);
321  Arc a3 = neg_gr.addArc(n2, n4);
322  Arc a4 = neg_gr.addArc(n3, n4);
323  Arc a5 = neg_gr.addArc(n3, n2);
324  Arc a6 = neg_gr.addArc(n5, n3);
325  Arc a7 = neg_gr.addArc(n5, n6);
326  Arc a8 = neg_gr.addArc(n6, n7);
327  Arc a9 = neg_gr.addArc(n7, n5);
328 
329  Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
330  ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
331  Digraph::NodeMap<int> neg_s(neg_gr, 0);
332 
333  neg_l2[a7] =  1000;
334  neg_l2[a8] = -1000;
335 
336  neg_s[n1] =  100;
337  neg_s[n4] = -100;
338 
339  neg_c[a1] =  100;
340  neg_c[a2] =   30;
341  neg_c[a3] =   20;
342  neg_c[a4] =   80;
343  neg_c[a5] =   50;
344  neg_c[a6] =   10;
345  neg_c[a7] =   80;
346  neg_c[a8] =   30;
347  neg_c[a9] = -120;
348
349  Digraph negs_gr;
350  Digraph::NodeMap<int> negs_s(negs_gr);
351  Digraph::ArcMap<int> negs_c(negs_gr);
352  ConstMap<Arc, int> negs_l(0), negs_u(1000);
353  n1 = negs_gr.addNode();
354  n2 = negs_gr.addNode();
355  negs_s[n1] = 100;
356  negs_s[n2] = -300;
357  negs_c[negs_gr.addArc(n1, n2)] = -1;
358
359
360  // A. Test NetworkSimplex with the default pivot rule
361  {
362    NetworkSimplex<Digraph> mcf(gr);
363
364    // Check the equality form
365    mcf.upperMap(u).costMap(c);
366    checkMcf(mcf, mcf.supplyMap(s1).run(),
367             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
368    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
369             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
370    mcf.lowerMap(l2);
371    checkMcf(mcf, mcf.supplyMap(s1).run(),
372             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
373    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
374             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
375    mcf.reset();
376    checkMcf(mcf, mcf.supplyMap(s1).run(),
377             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
378    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
379             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
380    mcf.reset();
381    checkMcf(mcf, mcf.run(),
382             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
383    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
384             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
385    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
386    checkMcf(mcf, mcf.run(),
387             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
388
389    // Check the GEQ form
390    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
391    checkMcf(mcf, mcf.run(),
392             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
393    mcf.supplyType(mcf.GEQ);
394    checkMcf(mcf, mcf.lowerMap(l2).run(),
395             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
396    mcf.supplyMap(s6);
397    checkMcf(mcf, mcf.run(),
398             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
399
400    // Check the LEQ form
401    mcf.reset().supplyType(mcf.LEQ);
402    mcf.upperMap(u).costMap(c).supplyMap(s6);
403    checkMcf(mcf, mcf.run(),
404             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
405    checkMcf(mcf, mcf.lowerMap(l2).run(),
406             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
407    mcf.supplyMap(s5);
408    checkMcf(mcf, mcf.run(),
409             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
410
411    // Check negative costs
412    NetworkSimplex<Digraph> neg_mcf(neg_gr);
413    neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
414    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
415      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
416    neg_mcf.upperMap(neg_u2);
417    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
418      neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
419    neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
420    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
421      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");
422     
423    NetworkSimplex<Digraph> negs_mcf(negs_gr);
424    negs_mcf.costMap(negs_c).supplyMap(negs_s);
425    checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
426      negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);
427  }
428
429  // B. Test NetworkSimplex with each pivot rule
430  {
431    NetworkSimplex<Digraph> mcf(gr);
432    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
433
434    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
435             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
436    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
437             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
438    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
439             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
440    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
441             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
442    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
443             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
444  }
445
446  return 0;
447}
Note: See TracBrowser for help on using the repository browser.