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 @ 998:7fdaa05a69a1

Last change on this file since 998:7fdaa05a69a1 was 877:141f9c0db4a3, checked in by Alpar Juttner <alpar@…>, 15 years ago

Unify the sources (#339)
File size: 16.8 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-2010
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
400template < typename MCF, typename Param >
401void runMcfLeqTests( Param param,
402                     const std::string &test_str = "" )
403{
404  // Tests for the LEQ form
405  MCF mcf1(gr);
406  mcf1.supplyType(mcf1.LEQ);
407  mcf1.upperMap(u).costMap(c).supplyMap(s6);
408  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6,
409           mcf1.OPTIMAL, true,   5080, test_str + "-19", LEQ);
410  mcf1.lowerMap(l2);
411  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
412           mcf1.OPTIMAL, true,   5930, test_str + "-20", LEQ);
413  mcf1.supplyMap(s5);
414  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
415           mcf1.INFEASIBLE, false,  0, test_str + "-21", LEQ);
416}
417
418
419int main()
420{
421  // Read the test networks
422  std::istringstream input(test_lgf);
423  DigraphReader<Digraph>(gr, input)
424    .arcMap("cost", c)
425    .arcMap("cap", u)
426    .arcMap("low1", l1)
427    .arcMap("low2", l2)
428    .arcMap("low3", l3)
429    .nodeMap("sup1", s1)
430    .nodeMap("sup2", s2)
431    .nodeMap("sup3", s3)
432    .nodeMap("sup4", s4)
433    .nodeMap("sup5", s5)
434    .nodeMap("sup6", s6)
435    .node("source", v)
436    .node("target", w)
437    .run();
438
439  std::istringstream neg_inp1(test_neg1_lgf);
440  DigraphReader<Digraph>(neg1_gr, neg_inp1)
441    .arcMap("cost", neg1_c)
442    .arcMap("low1", neg1_l1)
443    .arcMap("low2", neg1_l2)
444    .nodeMap("sup", neg1_s)
445    .run();
446
447  std::istringstream neg_inp2(test_neg2_lgf);
448  DigraphReader<Digraph>(neg2_gr, neg_inp2)
449    .arcMap("cost", neg2_c)
450    .nodeMap("sup", neg2_s)
451    .run();
452
453  // Check the interface of NetworkSimplex
454  {
455    typedef concepts::Digraph GR;
456    checkConcept< McfClassConcept<GR, int, int>,
457                  NetworkSimplex<GR> >();
458    checkConcept< McfClassConcept<GR, double, double>,
459                  NetworkSimplex<GR, double> >();
460    checkConcept< McfClassConcept<GR, int, double>,
461                  NetworkSimplex<GR, int, double> >();
462  }
463
464  // Check the interface of CapacityScaling
465  {
466    typedef concepts::Digraph GR;
467    checkConcept< McfClassConcept<GR, int, int>,
468                  CapacityScaling<GR> >();
469    checkConcept< McfClassConcept<GR, double, double>,
470                  CapacityScaling<GR, double> >();
471    checkConcept< McfClassConcept<GR, int, double>,
472                  CapacityScaling<GR, int, double> >();
473    typedef CapacityScaling<GR>::
474      SetHeap<concepts::Heap<int, RangeMap<int> > >::Create CAS;
475    checkConcept< McfClassConcept<GR, int, int>, CAS >();
476  }
477
478  // Check the interface of CostScaling
479  {
480    typedef concepts::Digraph GR;
481    checkConcept< McfClassConcept<GR, int, int>,
482                  CostScaling<GR> >();
483    checkConcept< McfClassConcept<GR, double, double>,
484                  CostScaling<GR, double> >();
485    checkConcept< McfClassConcept<GR, int, double>,
486                  CostScaling<GR, int, double> >();
487    typedef CostScaling<GR>::
488      SetLargeCost<double>::Create COS;
489    checkConcept< McfClassConcept<GR, int, int>, COS >();
490  }
491
492  // Check the interface of CycleCanceling
493  {
494    typedef concepts::Digraph GR;
495    checkConcept< McfClassConcept<GR, int, int>,
496                  CycleCanceling<GR> >();
497    checkConcept< McfClassConcept<GR, double, double>,
498                  CycleCanceling<GR, double> >();
499    checkConcept< McfClassConcept<GR, int, double>,
500                  CycleCanceling<GR, int, double> >();
501  }
502
503  // Test NetworkSimplex
504  {
505    typedef NetworkSimplex<Digraph> MCF;
506    runMcfGeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE", true);
507    runMcfLeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE");
508    runMcfGeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE", true);
509    runMcfLeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE");
510    runMcfGeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS", true);
511    runMcfLeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS");
512    runMcfGeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL", true);
513    runMcfLeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL");
514    runMcfGeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL", true);
515    runMcfLeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL");
516  }
517
518  // Test CapacityScaling
519  {
520    typedef CapacityScaling<Digraph> MCF;
521    runMcfGeqTests<MCF>(0, "SSP");
522    runMcfGeqTests<MCF>(2, "CAS");
523  }
524
525  // Test CostScaling
526  {
527    typedef CostScaling<Digraph> MCF;
528    runMcfGeqTests<MCF>(MCF::PUSH, "COS-PR");
529    runMcfGeqTests<MCF>(MCF::AUGMENT, "COS-AR");
530    runMcfGeqTests<MCF>(MCF::PARTIAL_AUGMENT, "COS-PAR");
531  }
532
533  // Test CycleCanceling
534  {
535    typedef CycleCanceling<Digraph> MCF;
536    runMcfGeqTests<MCF>(MCF::SIMPLE_CYCLE_CANCELING, "SCC");
537    runMcfGeqTests<MCF>(MCF::MINIMUM_MEAN_CYCLE_CANCELING, "MMCC");
538    runMcfGeqTests<MCF>(MCF::CANCEL_AND_TIGHTEN, "CAT");
539  }
540
541  return 0;
542}
Note: See TracBrowser for help on using the repository browser.