COIN-OR::LEMON - Graph Library

source: lemon/test/min_cost_flow_test.cc @ 1126:a30455cd0107

1.1 r1.1.4
Last change on this file since 1126:a30455cd0107 was 1081:f1398882a928, checked in by Alpar Juttner <alpar@…>, 13 years ago

Unify sources

File size: 13.9 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-2011
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      const Constraints& me = *this;
98
99      MCF mcf(me.g);
100      const MCF& const_mcf = mcf;
101
102      b = mcf.reset()
103             .lowerMap(me.lower)
104             .upperMap(me.upper)
105             .costMap(me.cost)
106             .supplyMap(me.sup)
107             .stSupply(me.n, me.n, me.k)
108             .run();
109
110      c = const_mcf.totalCost();
111      x = const_mcf.template totalCost<double>();
112      v = const_mcf.flow(me.a);
113      c = const_mcf.potential(me.n);
114      const_mcf.flowMap(fm);
115      const_mcf.potentialMap(pm);
116    }
117
118    typedef typename GR::Node Node;
119    typedef typename GR::Arc Arc;
120    typedef concepts::ReadMap<Node, Value> NM;
121    typedef concepts::ReadMap<Arc, Value> VAM;
122    typedef concepts::ReadMap<Arc, Cost> CAM;
123    typedef concepts::WriteMap<Arc, Value> FlowMap;
124    typedef concepts::WriteMap<Node, Cost> PotMap;
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;
134
135    FlowMap fm;
136    PotMap pm;
137    bool b;
138    double x;
139    typename MCF::Value v;
140    typename MCF::Cost c;
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,
150                const SM& supply, const FM& flow,
151                SupplyType type = EQ )
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];
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;
169  }
170
171  return true;
172}
173
174// Check the feasibility of the given potentials (dual soluiton)
175// using the "Complementary Slackness" optimality condition
176template < typename GR, typename LM, typename UM,
177           typename CM, typename SM, typename FM, typename PM >
178bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
179                     const CM& cost, const SM& supply, const FM& flow,
180                     const PM& pi, SupplyType type )
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  }
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];
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    }
204  }
205
206  return opt;
207}
208
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
243// Run a minimum cost flow algorithm and check the results
244template < typename MCF, typename GR,
245           typename LM, typename UM,
246           typename CM, typename SM,
247           typename PT >
248void checkMcf( const MCF& mcf, PT mcf_result,
249               const GR& gr, const LM& lower, const UM& upper,
250               const CM& cost, const SM& supply,
251               PT result, bool optimal, typename CM::Value total,
252               const std::string &test_id = "",
253               SupplyType type = EQ )
254{
255  check(mcf_result == result, "Wrong result " + test_id);
256  if (optimal) {
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),
262          "The flow is not feasible " + test_id);
263    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
264    check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
265          "Wrong potentials " + test_id);
266    check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
267          "Wrong dual cost " + test_id);
268  }
269}
270
271int main()
272{
273  // Check the interfaces
274  {
275    typedef concepts::Digraph GR;
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> >();
282  }
283
284  // Run various MCF tests
285  typedef ListDigraph Digraph;
286  DIGRAPH_TYPEDEFS(ListDigraph);
287
288  // Read the test digraph
289  Digraph gr;
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);
292  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
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)
301    .arcMap("low3", l3)
302    .nodeMap("sup1", s1)
303    .nodeMap("sup2", s2)
304    .nodeMap("sup3", s3)
305    .nodeMap("sup4", s4)
306    .nodeMap("sup5", s5)
307    .nodeMap("sup6", s6)
308    .node("source", v)
309    .node("target", w)
310    .run();
311
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();
321
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);
331
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);
335
336  neg_l2[a7] =  1000;
337  neg_l2[a8] = -1000;
338
339  neg_s[n1] =  100;
340  neg_s[n4] = -100;
341
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
362
363  // A. Test NetworkSimplex with the default pivot rule
364  {
365    NetworkSimplex<Digraph> mcf(gr);
366
367    // Check the equality form
368    mcf.upperMap(u).costMap(c);
369    checkMcf(mcf, mcf.supplyMap(s1).run(),
370             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
371    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
372             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
373    mcf.lowerMap(l2);
374    checkMcf(mcf, mcf.supplyMap(s1).run(),
375             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
376    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
377             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
378    mcf.reset();
379    checkMcf(mcf, mcf.supplyMap(s1).run(),
380             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
381    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
382             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
383    mcf.reset();
384    checkMcf(mcf, mcf.run(),
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");
391
392    // Check the GEQ form
393    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
394    checkMcf(mcf, mcf.run(),
395             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
396    mcf.supplyType(mcf.GEQ);
397    checkMcf(mcf, mcf.lowerMap(l2).run(),
398             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
399    mcf.supplyMap(s6);
400    checkMcf(mcf, mcf.run(),
401             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
402
403    // Check the LEQ form
404    mcf.reset().supplyType(mcf.LEQ);
405    mcf.upperMap(u).costMap(c).supplyMap(s6);
406    checkMcf(mcf, mcf.run(),
407             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
408    checkMcf(mcf, mcf.lowerMap(l2).run(),
409             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
410    mcf.supplyMap(s5);
411    checkMcf(mcf, mcf.run(),
412             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
413
414    // Check negative costs
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);
430  }
431
432  // B. Test NetworkSimplex with each pivot rule
433  {
434    NetworkSimplex<Digraph> mcf(gr);
435    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
436
437    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
438             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
439    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
440             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
441    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
442             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
443    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
444             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
445    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
446             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
447  }
448
449  return 0;
450}
Note: See TracBrowser for help on using the repository browser.