COIN-OR::LEMON - Graph Library

About

Downloads

Documentation

Development

Issue Tracker

Browse Tickets

Custom Query

Roadmap

Browse Source

Local Info

Mailing Lists

Search

source: lemon-main/test/min_cost_flow_test.cc @ 1160:bc571f16e1e9

Last change on this file since 1160:bc571f16e1e9 was 1118:ce1533650f7d, checked in by Alpar Juttner <alpar@…>, 11 years ago

Merge bugfix #474
File size: 17.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-2013
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#include <lemon/capacity_scaling.h>
28#include <lemon/cost_scaling.h>
29#include <lemon/cycle_canceling.h>
30
31#include <lemon/concepts/digraph.h>
32#include <lemon/concepts/heap.h>
33#include <lemon/concept_check.h>
34
35#include "test_tools.h"
36
37using namespace lemon;
38
39// Test networks
40char test_lgf[] =
41  "@nodes\n"
42  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
43  "    1    20   27    0   30   20   30\n"
44  "    2    -4    0    0    0   -8   -3\n"
45  "    3     0    0    0    0    0    0\n"
46  "    4     0    0    0    0    0    0\n"
47  "    5     9    0    0    0    6   11\n"
48  "    6    -6    0    0    0   -5   -6\n"
49  "    7     0    0    0    0    0    0\n"
50  "    8     0    0    0    0    0    3\n"
51  "    9     3    0    0    0    0    0\n"
52  "   10    -2    0    0    0   -7   -2\n"
53  "   11     0    0    0    0  -10    0\n"
54  "   12   -20  -27    0  -30  -30  -20\n"
55  "\n"
56  "@arcs\n"
57  "       cost  cap low1 low2 low3\n"
58  " 1  2    70   11    0    8    8\n"
59  " 1  3   150    3    0    1    0\n"
60  " 1  4    80   15    0    2    2\n"
61  " 2  8    80   12    0    0    0\n"
62  " 3  5   140    5    0    3    1\n"
63  " 4  6    60   10    0    1    0\n"
64  " 4  7    80    2    0    0    0\n"
65  " 4  8   110    3    0    0    0\n"
66  " 5  7    60   14    0    0    0\n"
67  " 5 11   120   12    0    0    0\n"
68  " 6  3     0    3    0    0    0\n"
69  " 6  9   140    4    0    0    0\n"
70  " 6 10    90    8    0    0    0\n"
71  " 7  1    30    5    0    0   -5\n"
72  " 8 12    60   16    0    4    3\n"
73  " 9 12    50    6    0    0    0\n"
74  "10 12    70   13    0    5    2\n"
75  "10  2   100    7    0    0    0\n"
76  "10  7    60   10    0    0   -3\n"
77  "11 10    20   14    0    6  -20\n"
78  "12 11    30   10    0    0  -10\n"
79  "\n"
80  "@attributes\n"
81  "source 1\n"
82  "target 12\n";
83
84char test_neg1_lgf[] =
85  "@nodes\n"
86  "label   sup\n"
87  "    1   100\n"
88  "    2     0\n"
89  "    3     0\n"
90  "    4  -100\n"
91  "    5     0\n"
92  "    6     0\n"
93  "    7     0\n"
94  "@arcs\n"
95  "      cost   low1   low2\n"
96  "1 2    100      0      0\n"
97  "1 3     30      0      0\n"
98  "2 4     20      0      0\n"
99  "3 4     80      0      0\n"
100  "3 2     50      0      0\n"
101  "5 3     10      0      0\n"
102  "5 6     80      0   1000\n"
103  "6 7     30      0  -1000\n"
104  "7 5   -120      0      0\n";
105
106char test_neg2_lgf[] =
107  "@nodes\n"
108  "label   sup\n"
109  "    1   100\n"
110  "    2  -300\n"
111  "@arcs\n"
112  "      cost\n"
113  "1 2     -1\n";
114
115
116// Test data
117typedef ListDigraph Digraph;
118DIGRAPH_TYPEDEFS(ListDigraph);
119
120Digraph gr;
121Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
122Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
123ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
124Node v, w;
125
126Digraph neg1_gr;
127Digraph::ArcMap<int> neg1_c(neg1_gr), neg1_l1(neg1_gr), neg1_l2(neg1_gr);
128ConstMap<Arc, int> neg1_u1(std::numeric_limits<int>::max()), neg1_u2(5000);
129Digraph::NodeMap<int> neg1_s(neg1_gr);
130
131Digraph neg2_gr;
132Digraph::ArcMap<int> neg2_c(neg2_gr);
133ConstMap<Arc, int> neg2_l(0), neg2_u(1000);
134Digraph::NodeMap<int> neg2_s(neg2_gr);
135
136
137enum SupplyType {
138  EQ,
139  GEQ,
140  LEQ
141};
142
143
144// Check the interface of an MCF algorithm
145template <typename GR, typename Value, typename Cost>
146class McfClassConcept
147{
148public:
149
150  template <typename MCF>
151  struct Constraints {
152    void constraints() {
153      checkConcept<concepts::Digraph, GR>();
154
155      const Constraints& me = *this;
156
157      MCF mcf(me.g);
158      const MCF& const_mcf = mcf;
159
160      b = mcf.reset().resetParams()
161             .lowerMap(me.lower)
162             .upperMap(me.upper)
163             .costMap(me.cost)
164             .supplyMap(me.sup)
165             .stSupply(me.n, me.n, me.k)
166             .run();
167
168      c = const_mcf.totalCost();
169      x = const_mcf.template totalCost<double>();
170      v = const_mcf.flow(me.a);
171      c = const_mcf.potential(me.n);
172      const_mcf.flowMap(fm);
173      const_mcf.potentialMap(pm);
174    }
175
176    typedef typename GR::Node Node;
177    typedef typename GR::Arc Arc;
178    typedef concepts::ReadMap<Node, Value> NM;
179    typedef concepts::ReadMap<Arc, Value> VAM;
180    typedef concepts::ReadMap<Arc, Cost> CAM;
181    typedef concepts::WriteMap<Arc, Value> FlowMap;
182    typedef concepts::WriteMap<Node, Cost> PotMap;
183
184    GR g;
185    VAM lower;
186    VAM upper;
187    CAM cost;
188    NM sup;
189    Node n;
190    Arc a;
191    Value k;
192
193    FlowMap fm;
194    PotMap pm;
195    bool b;
196    double x;
197    typename MCF::Value v;
198    typename MCF::Cost c;
199  };
200
201};
202
203
204// Check the feasibility of the given flow (primal soluiton)
205template < typename GR, typename LM, typename UM,
206           typename SM, typename FM >
207bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
208                const SM& supply, const FM& flow,
209                SupplyType type = EQ )
210{
211  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
212
213  for (ArcIt e(gr); e != INVALID; ++e) {
214    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
215  }
216
217  for (NodeIt n(gr); n != INVALID; ++n) {
218    typename SM::Value sum = 0;
219    for (OutArcIt e(gr, n); e != INVALID; ++e)
220      sum += flow[e];
221    for (InArcIt e(gr, n); e != INVALID; ++e)
222      sum -= flow[e];
223    bool b = (type ==  EQ && sum == supply[n]) ||
224             (type == GEQ && sum >= supply[n]) ||
225             (type == LEQ && sum <= supply[n]);
226    if (!b) return false;
227  }
228
229  return true;
230}
231
232// Check the feasibility of the given potentials (dual soluiton)
233// using the "Complementary Slackness" optimality condition
234template < typename GR, typename LM, typename UM,
235           typename CM, typename SM, typename FM, typename PM >
236bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
237                     const CM& cost, const SM& supply, const FM& flow,
238                     const PM& pi, SupplyType type )
239{
240  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
241
242  bool opt = true;
243  for (ArcIt e(gr); opt && e != INVALID; ++e) {
244    typename CM::Value red_cost =
245      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
246    opt = red_cost == 0 ||
247          (red_cost > 0 && flow[e] == lower[e]) ||
248          (red_cost < 0 && flow[e] == upper[e]);
249  }
250
251  for (NodeIt n(gr); opt && n != INVALID; ++n) {
252    typename SM::Value sum = 0;
253    for (OutArcIt e(gr, n); e != INVALID; ++e)
254      sum += flow[e];
255    for (InArcIt e(gr, n); e != INVALID; ++e)
256      sum -= flow[e];
257    if (type != LEQ) {
258      opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
259    } else {
260      opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
261    }
262  }
263
264  return opt;
265}
266
267// Check whether the dual cost is equal to the primal cost
268template < typename GR, typename LM, typename UM,
269           typename CM, typename SM, typename PM >
270bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
271                    const CM& cost, const SM& supply, const PM& pi,
272                    typename CM::Value total )
273{
274  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
275
276  typename CM::Value dual_cost = 0;
277  SM red_supply(gr);
278  for (NodeIt n(gr); n != INVALID; ++n) {
279    red_supply[n] = supply[n];
280  }
281  for (ArcIt a(gr); a != INVALID; ++a) {
282    if (lower[a] != 0) {
283      dual_cost += lower[a] * cost[a];
284      red_supply[gr.source(a)] -= lower[a];
285      red_supply[gr.target(a)] += lower[a];
286    }
287  }
288
289  for (NodeIt n(gr); n != INVALID; ++n) {
290    dual_cost -= red_supply[n] * pi[n];
291  }
292  for (ArcIt a(gr); a != INVALID; ++a) {
293    typename CM::Value red_cost =
294      cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
295    dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
296  }
297
298  return dual_cost == total;
299}
300
301// Run a minimum cost flow algorithm and check the results
302template < typename MCF, typename GR,
303           typename LM, typename UM,
304           typename CM, typename SM,
305           typename PT >
306void checkMcf( const MCF& mcf, PT mcf_result,
307               const GR& gr, const LM& lower, const UM& upper,
308               const CM& cost, const SM& supply,
309               PT result, bool optimal, typename CM::Value total,
310               const std::string &test_id = "",
311               SupplyType type = EQ )
312{
313  check(mcf_result == result, "Wrong result " + test_id);
314  if (optimal) {
315    typename GR::template ArcMap<typename SM::Value> flow(gr);
316    typename GR::template NodeMap<typename CM::Value> pi(gr);
317    mcf.flowMap(flow);
318    mcf.potentialMap(pi);
319    check(checkFlow(gr, lower, upper, supply, flow, type),
320          "The flow is not feasible " + test_id);
321    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
322    check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
323          "Wrong potentials " + test_id);
324    check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
325          "Wrong dual cost " + test_id);
326  }
327}
328
329template < typename MCF, typename Param >
330void runMcfGeqTests( Param param,
331                     const std::string &test_str = "",
332                     bool full_neg_cost_support = false )
333{
334  MCF mcf1(gr), mcf2(neg1_gr), mcf3(neg2_gr);
335
336  // Basic tests
337  mcf1.upperMap(u).costMap(c).supplyMap(s1);
338  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s1,
339           mcf1.OPTIMAL, true,     5240, test_str + "-1");
340  mcf1.stSupply(v, w, 27);
341  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s2,
342           mcf1.OPTIMAL, true,     7620, test_str + "-2");
343  mcf1.lowerMap(l2).supplyMap(s1);
344  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s1,
345           mcf1.OPTIMAL, true,     5970, test_str + "-3");
346  mcf1.stSupply(v, w, 27);
347  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s2,
348           mcf1.OPTIMAL, true,     8010, test_str + "-4");
349  mcf1.resetParams().supplyMap(s1);
350  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s1,
351           mcf1.OPTIMAL, true,       74, test_str + "-5");
352  mcf1.lowerMap(l2).stSupply(v, w, 27);
353  checkMcf(mcf1, mcf1.run(param), gr, l2, cu, cc, s2,
354           mcf1.OPTIMAL, true,       94, test_str + "-6");
355  mcf1.reset();
356  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s3,
357           mcf1.OPTIMAL, true,        0, test_str + "-7");
358  mcf1.lowerMap(l2).upperMap(u);
359  checkMcf(mcf1, mcf1.run(param), gr, l2, u, cc, s3,
360           mcf1.INFEASIBLE, false,    0, test_str + "-8");
361  mcf1.lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
362  checkMcf(mcf1, mcf1.run(param), gr, l3, u, c, s4,
363           mcf1.OPTIMAL, true,     6360, test_str + "-9");
364
365  // Tests for the GEQ form
366  mcf1.resetParams().upperMap(u).costMap(c).supplyMap(s5);
367  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s5,
368           mcf1.OPTIMAL, true,     3530, test_str + "-10", GEQ);
369  mcf1.lowerMap(l2);
370  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
371           mcf1.OPTIMAL, true,     4540, test_str + "-11", GEQ);
372  mcf1.supplyMap(s6);
373  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
374           mcf1.INFEASIBLE, false,    0, test_str + "-12", GEQ);
375
376  // Tests with negative costs
377  mcf2.lowerMap(neg1_l1).costMap(neg1_c).supplyMap(neg1_s);
378  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u1, neg1_c, neg1_s,
379           mcf2.UNBOUNDED, false,     0, test_str + "-13");
380  mcf2.upperMap(neg1_u2);
381  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u2, neg1_c, neg1_s,
382           mcf2.OPTIMAL, true,   -40000, test_str + "-14");
383  mcf2.resetParams().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s);
384  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l2, neg1_u1, neg1_c, neg1_s,
385           mcf2.UNBOUNDED, false,     0, test_str + "-15");
386
387  mcf3.costMap(neg2_c).supplyMap(neg2_s);
388  if (full_neg_cost_support) {
389    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
390             mcf3.OPTIMAL, true,   -300, test_str + "-16", GEQ);
391  } else {
392    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
393             mcf3.UNBOUNDED, false,   0, test_str + "-17", GEQ);
394  }
395  mcf3.upperMap(neg2_u);
396  checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
397           mcf3.OPTIMAL, true,     -300, test_str + "-18", GEQ);
398
399  // Tests for empty graph
400  Digraph gr0;
401  MCF mcf0(gr0);
402  mcf0.run(param);
403  check(mcf0.totalCost() == 0, "Wrong total cost"); 
404}
405
406template < typename MCF, typename Param >
407void runMcfLeqTests( Param param,
408                     const std::string &test_str = "" )
409{
410  // Tests for the LEQ form
411  MCF mcf1(gr);
412  mcf1.supplyType(mcf1.LEQ);
413  mcf1.upperMap(u).costMap(c).supplyMap(s6);
414  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6,
415           mcf1.OPTIMAL, true,   5080, test_str + "-19", LEQ);
416  mcf1.lowerMap(l2);
417  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
418           mcf1.OPTIMAL, true,   5930, test_str + "-20", LEQ);
419  mcf1.supplyMap(s5);
420  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
421           mcf1.INFEASIBLE, false,  0, test_str + "-21", LEQ);
422}
423
424
425int main()
426{
427  // Read the test networks
428  std::istringstream input(test_lgf);
429  DigraphReader<Digraph>(gr, input)
430    .arcMap("cost", c)
431    .arcMap("cap", u)
432    .arcMap("low1", l1)
433    .arcMap("low2", l2)
434    .arcMap("low3", l3)
435    .nodeMap("sup1", s1)
436    .nodeMap("sup2", s2)
437    .nodeMap("sup3", s3)
438    .nodeMap("sup4", s4)
439    .nodeMap("sup5", s5)
440    .nodeMap("sup6", s6)
441    .node("source", v)
442    .node("target", w)
443    .run();
444
445  std::istringstream neg_inp1(test_neg1_lgf);
446  DigraphReader<Digraph>(neg1_gr, neg_inp1)
447    .arcMap("cost", neg1_c)
448    .arcMap("low1", neg1_l1)
449    .arcMap("low2", neg1_l2)
450    .nodeMap("sup", neg1_s)
451    .run();
452
453  std::istringstream neg_inp2(test_neg2_lgf);
454  DigraphReader<Digraph>(neg2_gr, neg_inp2)
455    .arcMap("cost", neg2_c)
456    .nodeMap("sup", neg2_s)
457    .run();
458
459  // Check the interface of NetworkSimplex
460  {
461    typedef concepts::Digraph GR;
462    checkConcept< McfClassConcept<GR, int, int>,
463                  NetworkSimplex<GR> >();
464    checkConcept< McfClassConcept<GR, double, double>,
465                  NetworkSimplex<GR, double> >();
466    checkConcept< McfClassConcept<GR, int, double>,
467                  NetworkSimplex<GR, int, double> >();
468  }
469
470  // Check the interface of CapacityScaling
471  {
472    typedef concepts::Digraph GR;
473    checkConcept< McfClassConcept<GR, int, int>,
474                  CapacityScaling<GR> >();
475    checkConcept< McfClassConcept<GR, double, double>,
476                  CapacityScaling<GR, double> >();
477    checkConcept< McfClassConcept<GR, int, double>,
478                  CapacityScaling<GR, int, double> >();
479    typedef CapacityScaling<GR>::
480      SetHeap<concepts::Heap<int, RangeMap<int> > >::Create CAS;
481    checkConcept< McfClassConcept<GR, int, int>, CAS >();
482  }
483
484  // Check the interface of CostScaling
485  {
486    typedef concepts::Digraph GR;
487    checkConcept< McfClassConcept<GR, int, int>,
488                  CostScaling<GR> >();
489    checkConcept< McfClassConcept<GR, double, double>,
490                  CostScaling<GR, double> >();
491    checkConcept< McfClassConcept<GR, int, double>,
492                  CostScaling<GR, int, double> >();
493    typedef CostScaling<GR>::
494      SetLargeCost<double>::Create COS;
495    checkConcept< McfClassConcept<GR, int, int>, COS >();
496  }
497
498  // Check the interface of CycleCanceling
499  {
500    typedef concepts::Digraph GR;
501    checkConcept< McfClassConcept<GR, int, int>,
502                  CycleCanceling<GR> >();
503    checkConcept< McfClassConcept<GR, double, double>,
504                  CycleCanceling<GR, double> >();
505    checkConcept< McfClassConcept<GR, int, double>,
506                  CycleCanceling<GR, int, double> >();
507  }
508
509  // Test NetworkSimplex
510  {
511    typedef NetworkSimplex<Digraph> MCF;
512    runMcfGeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE", true);
513    runMcfLeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE");
514    runMcfGeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE", true);
515    runMcfLeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE");
516    runMcfGeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS", true);
517    runMcfLeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS");
518    runMcfGeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL", true);
519    runMcfLeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL");
520    runMcfGeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL", true);
521    runMcfLeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL");
522  }
523
524  // Test CapacityScaling
525  {
526    typedef CapacityScaling<Digraph> MCF;
527    runMcfGeqTests<MCF>(0, "SSP");
528    runMcfGeqTests<MCF>(2, "CAS");
529  }
530
531  // Test CostScaling
532  {
533    typedef CostScaling<Digraph> MCF;
534    runMcfGeqTests<MCF>(MCF::PUSH, "COS-PR");
535    runMcfGeqTests<MCF>(MCF::AUGMENT, "COS-AR");
536    runMcfGeqTests<MCF>(MCF::PARTIAL_AUGMENT, "COS-PAR");
537  }
538
539  // Test CycleCanceling
540  {
541    typedef CycleCanceling<Digraph> MCF;
542    runMcfGeqTests<MCF>(MCF::SIMPLE_CYCLE_CANCELING, "SCC");
543    runMcfGeqTests<MCF>(MCF::MINIMUM_MEAN_CYCLE_CANCELING, "MMCC");
544    runMcfGeqTests<MCF>(MCF::CANCEL_AND_TIGHTEN, "CAT");
545  }
546
547  return 0;
548}
Note: See TracBrowser for help on using the repository browser.