gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Merge and fix
0 8 2
merge default
8 files changed with 2191 insertions and 159 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
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
#ifndef LEMON_NETWORK_SIMPLEX_H
20
#define LEMON_NETWORK_SIMPLEX_H
21

	
22
/// \ingroup min_cost_flow
23
///
24
/// \file
25
/// \brief Network Simplex algorithm for finding a minimum cost flow.
26

	
27
#include <vector>
28
#include <limits>
29
#include <algorithm>
30

	
31
#include <lemon/core.h>
32
#include <lemon/math.h>
33
#include <lemon/maps.h>
34
#include <lemon/circulation.h>
35
#include <lemon/adaptors.h>
36

	
37
namespace lemon {
38

	
39
  /// \addtogroup min_cost_flow
40
  /// @{
41

	
42
  /// \brief Implementation of the primal Network Simplex algorithm
43
  /// for finding a \ref min_cost_flow "minimum cost flow".
44
  ///
45
  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
46
  /// for finding a \ref min_cost_flow "minimum cost flow".
47
  /// This algorithm is a specialized version of the linear programming
48
  /// simplex method directly for the minimum cost flow problem.
49
  /// It is one of the most efficient solution methods.
50
  ///
51
  /// In general this class is the fastest implementation available
52
  /// in LEMON for the minimum cost flow problem.
53
  /// Moreover it supports both direction of the supply/demand inequality
54
  /// constraints. For more information see \ref ProblemType.
55
  ///
56
  /// \tparam GR The digraph type the algorithm runs on.
57
  /// \tparam F The value type used for flow amounts, capacity bounds
58
  /// and supply values in the algorithm. By default it is \c int.
59
  /// \tparam C The value type used for costs and potentials in the
60
  /// algorithm. By default it is the same as \c F.
61
  ///
62
  /// \warning Both value types must be signed and all input data must
63
  /// be integer.
64
  ///
65
  /// \note %NetworkSimplex provides five different pivot rule
66
  /// implementations, from which the most efficient one is used
67
  /// by default. For more information see \ref PivotRule.
68
  template <typename GR, typename F = int, typename C = F>
69
  class NetworkSimplex
70
  {
71
  public:
72

	
73
    /// The flow type of the algorithm
74
    typedef F Flow;
75
    /// The cost type of the algorithm
76
    typedef C Cost;
77
#ifdef DOXYGEN
78
    /// The type of the flow map
79
    typedef GR::ArcMap<Flow> FlowMap;
80
    /// The type of the potential map
81
    typedef GR::NodeMap<Cost> PotentialMap;
82
#else
83
    /// The type of the flow map
84
    typedef typename GR::template ArcMap<Flow> FlowMap;
85
    /// The type of the potential map
86
    typedef typename GR::template NodeMap<Cost> PotentialMap;
87
#endif
88

	
89
  public:
90

	
91
    /// \brief Enum type for selecting the pivot rule.
92
    ///
93
    /// Enum type for selecting the pivot rule for the \ref run()
94
    /// function.
95
    ///
96
    /// \ref NetworkSimplex provides five different pivot rule
97
    /// implementations that significantly affect the running time
98
    /// of the algorithm.
99
    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
100
    /// proved to be the most efficient and the most robust on various
101
    /// test inputs according to our benchmark tests.
102
    /// However another pivot rule can be selected using the \ref run()
103
    /// function with the proper parameter.
104
    enum PivotRule {
105

	
106
      /// The First Eligible pivot rule.
107
      /// The next eligible arc is selected in a wraparound fashion
108
      /// in every iteration.
109
      FIRST_ELIGIBLE,
110

	
111
      /// The Best Eligible pivot rule.
112
      /// The best eligible arc is selected in every iteration.
113
      BEST_ELIGIBLE,
114

	
115
      /// The Block Search pivot rule.
116
      /// A specified number of arcs are examined in every iteration
117
      /// in a wraparound fashion and the best eligible arc is selected
118
      /// from this block.
119
      BLOCK_SEARCH,
120

	
121
      /// The Candidate List pivot rule.
122
      /// In a major iteration a candidate list is built from eligible arcs
123
      /// in a wraparound fashion and in the following minor iterations
124
      /// the best eligible arc is selected from this list.
125
      CANDIDATE_LIST,
126

	
127
      /// The Altering Candidate List pivot rule.
128
      /// It is a modified version of the Candidate List method.
129
      /// It keeps only the several best eligible arcs from the former
130
      /// candidate list and extends this list in every iteration.
131
      ALTERING_LIST
132
    };
133
    
134
    /// \brief Enum type for selecting the problem type.
135
    ///
136
    /// Enum type for selecting the problem type, i.e. the direction of
137
    /// the inequalities in the supply/demand constraints of the
138
    /// \ref min_cost_flow "minimum cost flow problem".
139
    ///
140
    /// The default problem type is \c GEQ, since this form is supported
141
    /// by other minimum cost flow algorithms and the \ref Circulation
142
    /// algorithm as well.
143
    /// The \c LEQ problem type can be selected using the \ref problemType()
144
    /// function.
145
    ///
146
    /// Note that the equality form is a special case of both problem type.
147
    enum ProblemType {
148

	
149
      /// This option means that there are "<em>greater or equal</em>"
150
      /// constraints in the defintion, i.e. the exact formulation of the
151
      /// problem is the following.
152
      /**
153
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
154
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
155
              sup(u) \quad \forall u\in V \f]
156
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
157
      */
158
      /// It means that the total demand must be greater or equal to the 
159
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
160
      /// negative) and all the supplies have to be carried out from 
161
      /// the supply nodes, but there could be demands that are not 
162
      /// satisfied.
163
      GEQ,
164
      /// It is just an alias for the \c GEQ option.
165
      CARRY_SUPPLIES = GEQ,
166

	
167
      /// This option means that there are "<em>less or equal</em>"
168
      /// constraints in the defintion, i.e. the exact formulation of the
169
      /// problem is the following.
170
      /**
171
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
172
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
173
              sup(u) \quad \forall u\in V \f]
174
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
175
      */
176
      /// It means that the total demand must be less or equal to the 
177
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
178
      /// positive) and all the demands have to be satisfied, but there
179
      /// could be supplies that are not carried out from the supply
180
      /// nodes.
181
      LEQ,
182
      /// It is just an alias for the \c LEQ option.
183
      SATISFY_DEMANDS = LEQ
184
    };
185

	
186
  private:
187

	
188
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
189

	
190
    typedef typename GR::template ArcMap<Flow> FlowArcMap;
191
    typedef typename GR::template ArcMap<Cost> CostArcMap;
192
    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
193

	
194
    typedef std::vector<Arc> ArcVector;
195
    typedef std::vector<Node> NodeVector;
196
    typedef std::vector<int> IntVector;
197
    typedef std::vector<bool> BoolVector;
198
    typedef std::vector<Flow> FlowVector;
199
    typedef std::vector<Cost> CostVector;
200

	
201
    // State constants for arcs
202
    enum ArcStateEnum {
203
      STATE_UPPER = -1,
204
      STATE_TREE  =  0,
205
      STATE_LOWER =  1
206
    };
207

	
208
  private:
209

	
210
    // Data related to the underlying digraph
211
    const GR &_graph;
212
    int _node_num;
213
    int _arc_num;
214

	
215
    // Parameters of the problem
216
    FlowArcMap *_plower;
217
    FlowArcMap *_pupper;
218
    CostArcMap *_pcost;
219
    FlowNodeMap *_psupply;
220
    bool _pstsup;
221
    Node _psource, _ptarget;
222
    Flow _pstflow;
223
    ProblemType _ptype;
224

	
225
    // Result maps
226
    FlowMap *_flow_map;
227
    PotentialMap *_potential_map;
228
    bool _local_flow;
229
    bool _local_potential;
230

	
231
    // Data structures for storing the digraph
232
    IntNodeMap _node_id;
233
    ArcVector _arc_ref;
234
    IntVector _source;
235
    IntVector _target;
236

	
237
    // Node and arc data
238
    FlowVector _cap;
239
    CostVector _cost;
240
    FlowVector _supply;
241
    FlowVector _flow;
242
    CostVector _pi;
243

	
244
    // Data for storing the spanning tree structure
245
    IntVector _parent;
246
    IntVector _pred;
247
    IntVector _thread;
248
    IntVector _rev_thread;
249
    IntVector _succ_num;
250
    IntVector _last_succ;
251
    IntVector _dirty_revs;
252
    BoolVector _forward;
253
    IntVector _state;
254
    int _root;
255

	
256
    // Temporary data used in the current pivot iteration
257
    int in_arc, join, u_in, v_in, u_out, v_out;
258
    int first, second, right, last;
259
    int stem, par_stem, new_stem;
260
    Flow delta;
261

	
262
  private:
263

	
264
    // Implementation of the First Eligible pivot rule
265
    class FirstEligiblePivotRule
266
    {
267
    private:
268

	
269
      // References to the NetworkSimplex class
270
      const IntVector  &_source;
271
      const IntVector  &_target;
272
      const CostVector &_cost;
273
      const IntVector  &_state;
274
      const CostVector &_pi;
275
      int &_in_arc;
276
      int _arc_num;
277

	
278
      // Pivot rule data
279
      int _next_arc;
280

	
281
    public:
282

	
283
      // Constructor
284
      FirstEligiblePivotRule(NetworkSimplex &ns) :
285
        _source(ns._source), _target(ns._target),
286
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
287
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
288
      {}
289

	
290
      // Find next entering arc
291
      bool findEnteringArc() {
292
        Cost c;
293
        for (int e = _next_arc; e < _arc_num; ++e) {
294
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
295
          if (c < 0) {
296
            _in_arc = e;
297
            _next_arc = e + 1;
298
            return true;
299
          }
300
        }
301
        for (int e = 0; e < _next_arc; ++e) {
302
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
303
          if (c < 0) {
304
            _in_arc = e;
305
            _next_arc = e + 1;
306
            return true;
307
          }
308
        }
309
        return false;
310
      }
311

	
312
    }; //class FirstEligiblePivotRule
313

	
314

	
315
    // Implementation of the Best Eligible pivot rule
316
    class BestEligiblePivotRule
317
    {
318
    private:
319

	
320
      // References to the NetworkSimplex class
321
      const IntVector  &_source;
322
      const IntVector  &_target;
323
      const CostVector &_cost;
324
      const IntVector  &_state;
325
      const CostVector &_pi;
326
      int &_in_arc;
327
      int _arc_num;
328

	
329
    public:
330

	
331
      // Constructor
332
      BestEligiblePivotRule(NetworkSimplex &ns) :
333
        _source(ns._source), _target(ns._target),
334
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
335
        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
336
      {}
337

	
338
      // Find next entering arc
339
      bool findEnteringArc() {
340
        Cost c, min = 0;
341
        for (int e = 0; e < _arc_num; ++e) {
342
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
343
          if (c < min) {
344
            min = c;
345
            _in_arc = e;
346
          }
347
        }
348
        return min < 0;
349
      }
350

	
351
    }; //class BestEligiblePivotRule
352

	
353

	
354
    // Implementation of the Block Search pivot rule
355
    class BlockSearchPivotRule
356
    {
357
    private:
358

	
359
      // References to the NetworkSimplex class
360
      const IntVector  &_source;
361
      const IntVector  &_target;
362
      const CostVector &_cost;
363
      const IntVector  &_state;
364
      const CostVector &_pi;
365
      int &_in_arc;
366
      int _arc_num;
367

	
368
      // Pivot rule data
369
      int _block_size;
370
      int _next_arc;
371

	
372
    public:
373

	
374
      // Constructor
375
      BlockSearchPivotRule(NetworkSimplex &ns) :
376
        _source(ns._source), _target(ns._target),
377
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
378
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
379
      {
380
        // The main parameters of the pivot rule
381
        const double BLOCK_SIZE_FACTOR = 2.0;
382
        const int MIN_BLOCK_SIZE = 10;
383

	
384
        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
385
                                MIN_BLOCK_SIZE );
386
      }
387

	
388
      // Find next entering arc
389
      bool findEnteringArc() {
390
        Cost c, min = 0;
391
        int cnt = _block_size;
392
        int e, min_arc = _next_arc;
393
        for (e = _next_arc; e < _arc_num; ++e) {
394
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
395
          if (c < min) {
396
            min = c;
397
            min_arc = e;
398
          }
399
          if (--cnt == 0) {
400
            if (min < 0) break;
401
            cnt = _block_size;
402
          }
403
        }
404
        if (min == 0 || cnt > 0) {
405
          for (e = 0; e < _next_arc; ++e) {
406
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
407
            if (c < min) {
408
              min = c;
409
              min_arc = e;
410
            }
411
            if (--cnt == 0) {
412
              if (min < 0) break;
413
              cnt = _block_size;
414
            }
415
          }
416
        }
417
        if (min >= 0) return false;
418
        _in_arc = min_arc;
419
        _next_arc = e;
420
        return true;
421
      }
422

	
423
    }; //class BlockSearchPivotRule
424

	
425

	
426
    // Implementation of the Candidate List pivot rule
427
    class CandidateListPivotRule
428
    {
429
    private:
430

	
431
      // References to the NetworkSimplex class
432
      const IntVector  &_source;
433
      const IntVector  &_target;
434
      const CostVector &_cost;
435
      const IntVector  &_state;
436
      const CostVector &_pi;
437
      int &_in_arc;
438
      int _arc_num;
439

	
440
      // Pivot rule data
441
      IntVector _candidates;
442
      int _list_length, _minor_limit;
443
      int _curr_length, _minor_count;
444
      int _next_arc;
445

	
446
    public:
447

	
448
      /// Constructor
449
      CandidateListPivotRule(NetworkSimplex &ns) :
450
        _source(ns._source), _target(ns._target),
451
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
452
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
453
      {
454
        // The main parameters of the pivot rule
455
        const double LIST_LENGTH_FACTOR = 1.0;
456
        const int MIN_LIST_LENGTH = 10;
457
        const double MINOR_LIMIT_FACTOR = 0.1;
458
        const int MIN_MINOR_LIMIT = 3;
459

	
460
        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
461
                                 MIN_LIST_LENGTH );
462
        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
463
                                 MIN_MINOR_LIMIT );
464
        _curr_length = _minor_count = 0;
465
        _candidates.resize(_list_length);
466
      }
467

	
468
      /// Find next entering arc
469
      bool findEnteringArc() {
470
        Cost min, c;
471
        int e, min_arc = _next_arc;
472
        if (_curr_length > 0 && _minor_count < _minor_limit) {
473
          // Minor iteration: select the best eligible arc from the
474
          // current candidate list
475
          ++_minor_count;
476
          min = 0;
477
          for (int i = 0; i < _curr_length; ++i) {
478
            e = _candidates[i];
479
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
480
            if (c < min) {
481
              min = c;
482
              min_arc = e;
483
            }
484
            if (c >= 0) {
485
              _candidates[i--] = _candidates[--_curr_length];
486
            }
487
          }
488
          if (min < 0) {
489
            _in_arc = min_arc;
490
            return true;
491
          }
492
        }
493

	
494
        // Major iteration: build a new candidate list
495
        min = 0;
496
        _curr_length = 0;
497
        for (e = _next_arc; e < _arc_num; ++e) {
498
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
499
          if (c < 0) {
500
            _candidates[_curr_length++] = e;
501
            if (c < min) {
502
              min = c;
503
              min_arc = e;
504
            }
505
            if (_curr_length == _list_length) break;
506
          }
507
        }
508
        if (_curr_length < _list_length) {
509
          for (e = 0; e < _next_arc; ++e) {
510
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
511
            if (c < 0) {
512
              _candidates[_curr_length++] = e;
513
              if (c < min) {
514
                min = c;
515
                min_arc = e;
516
              }
517
              if (_curr_length == _list_length) break;
518
            }
519
          }
520
        }
521
        if (_curr_length == 0) return false;
522
        _minor_count = 1;
523
        _in_arc = min_arc;
524
        _next_arc = e;
525
        return true;
526
      }
527

	
528
    }; //class CandidateListPivotRule
529

	
530

	
531
    // Implementation of the Altering Candidate List pivot rule
532
    class AlteringListPivotRule
533
    {
534
    private:
535

	
536
      // References to the NetworkSimplex class
537
      const IntVector  &_source;
538
      const IntVector  &_target;
539
      const CostVector &_cost;
540
      const IntVector  &_state;
541
      const CostVector &_pi;
542
      int &_in_arc;
543
      int _arc_num;
544

	
545
      // Pivot rule data
546
      int _block_size, _head_length, _curr_length;
547
      int _next_arc;
548
      IntVector _candidates;
549
      CostVector _cand_cost;
550

	
551
      // Functor class to compare arcs during sort of the candidate list
552
      class SortFunc
553
      {
554
      private:
555
        const CostVector &_map;
556
      public:
557
        SortFunc(const CostVector &map) : _map(map) {}
558
        bool operator()(int left, int right) {
559
          return _map[left] > _map[right];
560
        }
561
      };
562

	
563
      SortFunc _sort_func;
564

	
565
    public:
566

	
567
      // Constructor
568
      AlteringListPivotRule(NetworkSimplex &ns) :
569
        _source(ns._source), _target(ns._target),
570
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
571
        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
572
        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
573
      {
574
        // The main parameters of the pivot rule
575
        const double BLOCK_SIZE_FACTOR = 1.5;
576
        const int MIN_BLOCK_SIZE = 10;
577
        const double HEAD_LENGTH_FACTOR = 0.1;
578
        const int MIN_HEAD_LENGTH = 3;
579

	
580
        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
581
                                MIN_BLOCK_SIZE );
582
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
583
                                 MIN_HEAD_LENGTH );
584
        _candidates.resize(_head_length + _block_size);
585
        _curr_length = 0;
586
      }
587

	
588
      // Find next entering arc
589
      bool findEnteringArc() {
590
        // Check the current candidate list
591
        int e;
592
        for (int i = 0; i < _curr_length; ++i) {
593
          e = _candidates[i];
594
          _cand_cost[e] = _state[e] *
595
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
596
          if (_cand_cost[e] >= 0) {
597
            _candidates[i--] = _candidates[--_curr_length];
598
          }
599
        }
600

	
601
        // Extend the list
602
        int cnt = _block_size;
603
        int last_arc = 0;
604
        int limit = _head_length;
605

	
606
        for (int e = _next_arc; e < _arc_num; ++e) {
607
          _cand_cost[e] = _state[e] *
608
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
609
          if (_cand_cost[e] < 0) {
610
            _candidates[_curr_length++] = e;
611
            last_arc = e;
612
          }
613
          if (--cnt == 0) {
614
            if (_curr_length > limit) break;
615
            limit = 0;
616
            cnt = _block_size;
617
          }
618
        }
619
        if (_curr_length <= limit) {
620
          for (int e = 0; e < _next_arc; ++e) {
621
            _cand_cost[e] = _state[e] *
622
              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
623
            if (_cand_cost[e] < 0) {
624
              _candidates[_curr_length++] = e;
625
              last_arc = e;
626
            }
627
            if (--cnt == 0) {
628
              if (_curr_length > limit) break;
629
              limit = 0;
630
              cnt = _block_size;
631
            }
632
          }
633
        }
634
        if (_curr_length == 0) return false;
635
        _next_arc = last_arc + 1;
636

	
637
        // Make heap of the candidate list (approximating a partial sort)
638
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
639
                   _sort_func );
640

	
641
        // Pop the first element of the heap
642
        _in_arc = _candidates[0];
643
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
644
                  _sort_func );
645
        _curr_length = std::min(_head_length, _curr_length - 1);
646
        return true;
647
      }
648

	
649
    }; //class AlteringListPivotRule
650

	
651
  public:
652

	
653
    /// \brief Constructor.
654
    ///
655
    /// The constructor of the class.
656
    ///
657
    /// \param graph The digraph the algorithm runs on.
658
    NetworkSimplex(const GR& graph) :
659
      _graph(graph),
660
      _plower(NULL), _pupper(NULL), _pcost(NULL),
661
      _psupply(NULL), _pstsup(false), _ptype(GEQ),
662
      _flow_map(NULL), _potential_map(NULL),
663
      _local_flow(false), _local_potential(false),
664
      _node_id(graph)
665
    {
666
      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
667
                   std::numeric_limits<Flow>::is_signed,
668
        "The flow type of NetworkSimplex must be signed integer");
669
      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
670
                   std::numeric_limits<Cost>::is_signed,
671
        "The cost type of NetworkSimplex must be signed integer");
672
    }
673

	
674
    /// Destructor.
675
    ~NetworkSimplex() {
676
      if (_local_flow) delete _flow_map;
677
      if (_local_potential) delete _potential_map;
678
    }
679

	
680
    /// \name Parameters
681
    /// The parameters of the algorithm can be specified using these
682
    /// functions.
683

	
684
    /// @{
685

	
686
    /// \brief Set the lower bounds on the arcs.
687
    ///
688
    /// This function sets the lower bounds on the arcs.
689
    /// If neither this function nor \ref boundMaps() is used before
690
    /// calling \ref run(), the lower bounds will be set to zero
691
    /// on all arcs.
692
    ///
693
    /// \param map An arc map storing the lower bounds.
694
    /// Its \c Value type must be convertible to the \c Flow type
695
    /// of the algorithm.
696
    ///
697
    /// \return <tt>(*this)</tt>
698
    template <typename LOWER>
699
    NetworkSimplex& lowerMap(const LOWER& map) {
700
      delete _plower;
701
      _plower = new FlowArcMap(_graph);
702
      for (ArcIt a(_graph); a != INVALID; ++a) {
703
        (*_plower)[a] = map[a];
704
      }
705
      return *this;
706
    }
707

	
708
    /// \brief Set the upper bounds (capacities) on the arcs.
709
    ///
710
    /// This function sets the upper bounds (capacities) on the arcs.
711
    /// If none of the functions \ref upperMap(), \ref capacityMap()
712
    /// and \ref boundMaps() is used before calling \ref run(),
713
    /// the upper bounds (capacities) will be set to
714
    /// \c std::numeric_limits<Flow>::max() on all arcs.
715
    ///
716
    /// \param map An arc map storing the upper bounds.
717
    /// Its \c Value type must be convertible to the \c Flow type
718
    /// of the algorithm.
719
    ///
720
    /// \return <tt>(*this)</tt>
721
    template<typename UPPER>
722
    NetworkSimplex& upperMap(const UPPER& map) {
723
      delete _pupper;
724
      _pupper = new FlowArcMap(_graph);
725
      for (ArcIt a(_graph); a != INVALID; ++a) {
726
        (*_pupper)[a] = map[a];
727
      }
728
      return *this;
729
    }
730

	
731
    /// \brief Set the upper bounds (capacities) on the arcs.
732
    ///
733
    /// This function sets the upper bounds (capacities) on the arcs.
734
    /// It is just an alias for \ref upperMap().
735
    ///
736
    /// \return <tt>(*this)</tt>
737
    template<typename CAP>
738
    NetworkSimplex& capacityMap(const CAP& map) {
739
      return upperMap(map);
740
    }
741

	
742
    /// \brief Set the lower and upper bounds on the arcs.
743
    ///
744
    /// This function sets the lower and upper bounds on the arcs.
745
    /// If neither this function nor \ref lowerMap() is used before
746
    /// calling \ref run(), the lower bounds will be set to zero
747
    /// on all arcs.
748
    /// If none of the functions \ref upperMap(), \ref capacityMap()
749
    /// and \ref boundMaps() is used before calling \ref run(),
750
    /// the upper bounds (capacities) will be set to
751
    /// \c std::numeric_limits<Flow>::max() on all arcs.
752
    ///
753
    /// \param lower An arc map storing the lower bounds.
754
    /// \param upper An arc map storing the upper bounds.
755
    ///
756
    /// The \c Value type of the maps must be convertible to the
757
    /// \c Flow type of the algorithm.
758
    ///
759
    /// \note This function is just a shortcut of calling \ref lowerMap()
760
    /// and \ref upperMap() separately.
761
    ///
762
    /// \return <tt>(*this)</tt>
763
    template <typename LOWER, typename UPPER>
764
    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
765
      return lowerMap(lower).upperMap(upper);
766
    }
767

	
768
    /// \brief Set the costs of the arcs.
769
    ///
770
    /// This function sets the costs of the arcs.
771
    /// If it is not used before calling \ref run(), the costs
772
    /// will be set to \c 1 on all arcs.
773
    ///
774
    /// \param map An arc map storing the costs.
775
    /// Its \c Value type must be convertible to the \c Cost type
776
    /// of the algorithm.
777
    ///
778
    /// \return <tt>(*this)</tt>
779
    template<typename COST>
780
    NetworkSimplex& costMap(const COST& map) {
781
      delete _pcost;
782
      _pcost = new CostArcMap(_graph);
783
      for (ArcIt a(_graph); a != INVALID; ++a) {
784
        (*_pcost)[a] = map[a];
785
      }
786
      return *this;
787
    }
788

	
789
    /// \brief Set the supply values of the nodes.
790
    ///
791
    /// This function sets the supply values of the nodes.
792
    /// If neither this function nor \ref stSupply() is used before
793
    /// calling \ref run(), the supply of each node will be set to zero.
794
    /// (It makes sense only if non-zero lower bounds are given.)
795
    ///
796
    /// \param map A node map storing the supply values.
797
    /// Its \c Value type must be convertible to the \c Flow type
798
    /// of the algorithm.
799
    ///
800
    /// \return <tt>(*this)</tt>
801
    template<typename SUP>
802
    NetworkSimplex& supplyMap(const SUP& map) {
803
      delete _psupply;
804
      _pstsup = false;
805
      _psupply = new FlowNodeMap(_graph);
806
      for (NodeIt n(_graph); n != INVALID; ++n) {
807
        (*_psupply)[n] = map[n];
808
      }
809
      return *this;
810
    }
811

	
812
    /// \brief Set single source and target nodes and a supply value.
813
    ///
814
    /// This function sets a single source node and a single target node
815
    /// and the required flow value.
816
    /// If neither this function nor \ref supplyMap() is used before
817
    /// calling \ref run(), the supply of each node will be set to zero.
818
    /// (It makes sense only if non-zero lower bounds are given.)
819
    ///
820
    /// \param s The source node.
821
    /// \param t The target node.
822
    /// \param k The required amount of flow from node \c s to node \c t
823
    /// (i.e. the supply of \c s and the demand of \c t).
824
    ///
825
    /// \return <tt>(*this)</tt>
826
    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
827
      delete _psupply;
828
      _psupply = NULL;
829
      _pstsup = true;
830
      _psource = s;
831
      _ptarget = t;
832
      _pstflow = k;
833
      return *this;
834
    }
835
    
836
    /// \brief Set the problem type.
837
    ///
838
    /// This function sets the problem type for the algorithm.
839
    /// If it is not used before calling \ref run(), the \ref GEQ problem
840
    /// type will be used.
841
    ///
842
    /// For more information see \ref ProblemType.
843
    ///
844
    /// \return <tt>(*this)</tt>
845
    NetworkSimplex& problemType(ProblemType problem_type) {
846
      _ptype = problem_type;
847
      return *this;
848
    }
849

	
850
    /// \brief Set the flow map.
851
    ///
852
    /// This function sets the flow map.
853
    /// If it is not used before calling \ref run(), an instance will
854
    /// be allocated automatically. The destructor deallocates this
855
    /// automatically allocated map, of course.
856
    ///
857
    /// \return <tt>(*this)</tt>
858
    NetworkSimplex& flowMap(FlowMap& map) {
859
      if (_local_flow) {
860
        delete _flow_map;
861
        _local_flow = false;
862
      }
863
      _flow_map = &map;
864
      return *this;
865
    }
866

	
867
    /// \brief Set the potential map.
868
    ///
869
    /// This function sets the potential map, which is used for storing
870
    /// the dual solution.
871
    /// If it is not used before calling \ref run(), an instance will
872
    /// be allocated automatically. The destructor deallocates this
873
    /// automatically allocated map, of course.
874
    ///
875
    /// \return <tt>(*this)</tt>
876
    NetworkSimplex& potentialMap(PotentialMap& map) {
877
      if (_local_potential) {
878
        delete _potential_map;
879
        _local_potential = false;
880
      }
881
      _potential_map = &map;
882
      return *this;
883
    }
884
    
885
    /// @}
886

	
887
    /// \name Execution Control
888
    /// The algorithm can be executed using \ref run().
889

	
890
    /// @{
891

	
892
    /// \brief Run the algorithm.
893
    ///
894
    /// This function runs the algorithm.
895
    /// The paramters can be specified using functions \ref lowerMap(),
896
    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
897
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
898
    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
899
    /// For example,
900
    /// \code
901
    ///   NetworkSimplex<ListDigraph> ns(graph);
902
    ///   ns.boundMaps(lower, upper).costMap(cost)
903
    ///     .supplyMap(sup).run();
904
    /// \endcode
905
    ///
906
    /// This function can be called more than once. All the parameters
907
    /// that have been given are kept for the next call, unless
908
    /// \ref reset() is called, thus only the modified parameters
909
    /// have to be set again. See \ref reset() for examples.
910
    ///
911
    /// \param pivot_rule The pivot rule that will be used during the
912
    /// algorithm. For more information see \ref PivotRule.
913
    ///
914
    /// \return \c true if a feasible flow can be found.
915
    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
916
      return init() && start(pivot_rule);
917
    }
918

	
919
    /// \brief Reset all the parameters that have been given before.
920
    ///
921
    /// This function resets all the paramaters that have been given
922
    /// before using functions \ref lowerMap(), \ref upperMap(),
923
    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
924
    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
925
    /// \ref flowMap() and \ref potentialMap().
926
    ///
927
    /// It is useful for multiple run() calls. If this function is not
928
    /// used, all the parameters given before are kept for the next
929
    /// \ref run() call.
930
    ///
931
    /// For example,
932
    /// \code
933
    ///   NetworkSimplex<ListDigraph> ns(graph);
934
    ///
935
    ///   // First run
936
    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
937
    ///     .supplyMap(sup).run();
938
    ///
939
    ///   // Run again with modified cost map (reset() is not called,
940
    ///   // so only the cost map have to be set again)
941
    ///   cost[e] += 100;
942
    ///   ns.costMap(cost).run();
943
    ///
944
    ///   // Run again from scratch using reset()
945
    ///   // (the lower bounds will be set to zero on all arcs)
946
    ///   ns.reset();
947
    ///   ns.capacityMap(cap).costMap(cost)
948
    ///     .supplyMap(sup).run();
949
    /// \endcode
950
    ///
951
    /// \return <tt>(*this)</tt>
952
    NetworkSimplex& reset() {
953
      delete _plower;
954
      delete _pupper;
955
      delete _pcost;
956
      delete _psupply;
957
      _plower = NULL;
958
      _pupper = NULL;
959
      _pcost = NULL;
960
      _psupply = NULL;
961
      _pstsup = false;
962
      _ptype = GEQ;
963
      if (_local_flow) delete _flow_map;
964
      if (_local_potential) delete _potential_map;
965
      _flow_map = NULL;
966
      _potential_map = NULL;
967
      _local_flow = false;
968
      _local_potential = false;
969

	
970
      return *this;
971
    }
972

	
973
    /// @}
974

	
975
    /// \name Query Functions
976
    /// The results of the algorithm can be obtained using these
977
    /// functions.\n
978
    /// The \ref run() function must be called before using them.
979

	
980
    /// @{
981

	
982
    /// \brief Return the total cost of the found flow.
983
    ///
984
    /// This function returns the total cost of the found flow.
985
    /// The complexity of the function is O(e).
986
    ///
987
    /// \note The return type of the function can be specified as a
988
    /// template parameter. For example,
989
    /// \code
990
    ///   ns.totalCost<double>();
991
    /// \endcode
992
    /// It is useful if the total cost cannot be stored in the \c Cost
993
    /// type of the algorithm, which is the default return type of the
994
    /// function.
995
    ///
996
    /// \pre \ref run() must be called before using this function.
997
    template <typename Num>
998
    Num totalCost() const {
999
      Num c = 0;
1000
      if (_pcost) {
1001
        for (ArcIt e(_graph); e != INVALID; ++e)
1002
          c += (*_flow_map)[e] * (*_pcost)[e];
1003
      } else {
1004
        for (ArcIt e(_graph); e != INVALID; ++e)
1005
          c += (*_flow_map)[e];
1006
      }
1007
      return c;
1008
    }
1009

	
1010
#ifndef DOXYGEN
1011
    Cost totalCost() const {
1012
      return totalCost<Cost>();
1013
    }
1014
#endif
1015

	
1016
    /// \brief Return the flow on the given arc.
1017
    ///
1018
    /// This function returns the flow on the given arc.
1019
    ///
1020
    /// \pre \ref run() must be called before using this function.
1021
    Flow flow(const Arc& a) const {
1022
      return (*_flow_map)[a];
1023
    }
1024

	
1025
    /// \brief Return a const reference to the flow map.
1026
    ///
1027
    /// This function returns a const reference to an arc map storing
1028
    /// the found flow.
1029
    ///
1030
    /// \pre \ref run() must be called before using this function.
1031
    const FlowMap& flowMap() const {
1032
      return *_flow_map;
1033
    }
1034

	
1035
    /// \brief Return the potential (dual value) of the given node.
1036
    ///
1037
    /// This function returns the potential (dual value) of the
1038
    /// given node.
1039
    ///
1040
    /// \pre \ref run() must be called before using this function.
1041
    Cost potential(const Node& n) const {
1042
      return (*_potential_map)[n];
1043
    }
1044

	
1045
    /// \brief Return a const reference to the potential map
1046
    /// (the dual solution).
1047
    ///
1048
    /// This function returns a const reference to a node map storing
1049
    /// the found potentials, which form the dual solution of the
1050
    /// \ref min_cost_flow "minimum cost flow" problem.
1051
    ///
1052
    /// \pre \ref run() must be called before using this function.
1053
    const PotentialMap& potentialMap() const {
1054
      return *_potential_map;
1055
    }
1056

	
1057
    /// @}
1058

	
1059
  private:
1060

	
1061
    // Initialize internal data structures
1062
    bool init() {
1063
      // Initialize result maps
1064
      if (!_flow_map) {
1065
        _flow_map = new FlowMap(_graph);
1066
        _local_flow = true;
1067
      }
1068
      if (!_potential_map) {
1069
        _potential_map = new PotentialMap(_graph);
1070
        _local_potential = true;
1071
      }
1072

	
1073
      // Initialize vectors
1074
      _node_num = countNodes(_graph);
1075
      _arc_num = countArcs(_graph);
1076
      int all_node_num = _node_num + 1;
1077
      int all_arc_num = _arc_num + _node_num;
1078
      if (_node_num == 0) return false;
1079

	
1080
      _arc_ref.resize(_arc_num);
1081
      _source.resize(all_arc_num);
1082
      _target.resize(all_arc_num);
1083

	
1084
      _cap.resize(all_arc_num);
1085
      _cost.resize(all_arc_num);
1086
      _supply.resize(all_node_num);
1087
      _flow.resize(all_arc_num);
1088
      _pi.resize(all_node_num);
1089

	
1090
      _parent.resize(all_node_num);
1091
      _pred.resize(all_node_num);
1092
      _forward.resize(all_node_num);
1093
      _thread.resize(all_node_num);
1094
      _rev_thread.resize(all_node_num);
1095
      _succ_num.resize(all_node_num);
1096
      _last_succ.resize(all_node_num);
1097
      _state.resize(all_arc_num);
1098

	
1099
      // Initialize node related data
1100
      bool valid_supply = true;
1101
      Flow sum_supply = 0;
1102
      if (!_pstsup && !_psupply) {
1103
        _pstsup = true;
1104
        _psource = _ptarget = NodeIt(_graph);
1105
        _pstflow = 0;
1106
      }
1107
      if (_psupply) {
1108
        int i = 0;
1109
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1110
          _node_id[n] = i;
1111
          _supply[i] = (*_psupply)[n];
1112
          sum_supply += _supply[i];
1113
        }
1114
        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
1115
                       (_ptype == LEQ && sum_supply >= 0);
1116
      } else {
1117
        int i = 0;
1118
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1119
          _node_id[n] = i;
1120
          _supply[i] = 0;
1121
        }
1122
        _supply[_node_id[_psource]] =  _pstflow;
1123
        _supply[_node_id[_ptarget]] = -_pstflow;
1124
      }
1125
      if (!valid_supply) return false;
1126

	
1127
      // Infinite capacity value
1128
      Flow inf_cap =
1129
        std::numeric_limits<Flow>::has_infinity ?
1130
        std::numeric_limits<Flow>::infinity() :
1131
        std::numeric_limits<Flow>::max();
1132

	
1133
      // Initialize artifical cost
1134
      Cost art_cost;
1135
      if (std::numeric_limits<Cost>::is_exact) {
1136
        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1137
      } else {
1138
        art_cost = std::numeric_limits<Cost>::min();
1139
        for (int i = 0; i != _arc_num; ++i) {
1140
          if (_cost[i] > art_cost) art_cost = _cost[i];
1141
        }
1142
        art_cost = (art_cost + 1) * _node_num;
1143
      }
1144

	
1145
      // Run Circulation to check if a feasible solution exists
1146
      typedef ConstMap<Arc, Flow> ConstArcMap;
1147
      FlowNodeMap *csup = NULL;
1148
      bool local_csup = false;
1149
      if (_psupply) {
1150
        csup = _psupply;
1151
      } else {
1152
        csup = new FlowNodeMap(_graph, 0);
1153
        (*csup)[_psource] =  _pstflow;
1154
        (*csup)[_ptarget] = -_pstflow;
1155
        local_csup = true;
1156
      }
1157
      bool circ_result = false;
1158
      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
1159
        // GEQ problem type
1160
        if (_plower) {
1161
          if (_pupper) {
1162
            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
1163
              circ(_graph, *_plower, *_pupper, *csup);
1164
            circ_result = circ.run();
1165
          } else {
1166
            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
1167
              circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
1168
            circ_result = circ.run();
1169
          }
1170
        } else {
1171
          if (_pupper) {
1172
            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
1173
              circ(_graph, ConstArcMap(0), *_pupper, *csup);
1174
            circ_result = circ.run();
1175
          } else {
1176
            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
1177
              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
1178
            circ_result = circ.run();
1179
          }
1180
        }
1181
      } else {
1182
        // LEQ problem type
1183
        typedef ReverseDigraph<const GR> RevGraph;
1184
        typedef NegMap<FlowNodeMap> NegNodeMap;
1185
        RevGraph rgraph(_graph);
1186
        NegNodeMap neg_csup(*csup);
1187
        if (_plower) {
1188
          if (_pupper) {
1189
            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
1190
              circ(rgraph, *_plower, *_pupper, neg_csup);
1191
            circ_result = circ.run();
1192
          } else {
1193
            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
1194
              circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
1195
            circ_result = circ.run();
1196
          }
1197
        } else {
1198
          if (_pupper) {
1199
            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
1200
              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
1201
            circ_result = circ.run();
1202
          } else {
1203
            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
1204
              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
1205
            circ_result = circ.run();
1206
          }
1207
        }
1208
      }
1209
      if (local_csup) delete csup;
1210
      if (!circ_result) return false;
1211

	
1212
      // Set data for the artificial root node
1213
      _root = _node_num;
1214
      _parent[_root] = -1;
1215
      _pred[_root] = -1;
1216
      _thread[_root] = 0;
1217
      _rev_thread[0] = _root;
1218
      _succ_num[_root] = all_node_num;
1219
      _last_succ[_root] = _root - 1;
1220
      _supply[_root] = -sum_supply;
1221
      if (sum_supply < 0) {
1222
        _pi[_root] = -art_cost;
1223
      } else {
1224
        _pi[_root] = art_cost;
1225
      }
1226

	
1227
      // Store the arcs in a mixed order
1228
      int k = std::max(int(sqrt(_arc_num)), 10);
1229
      int i = 0;
1230
      for (ArcIt e(_graph); e != INVALID; ++e) {
1231
        _arc_ref[i] = e;
1232
        if ((i += k) >= _arc_num) i = (i % k) + 1;
1233
      }
1234

	
1235
      // Initialize arc maps
1236
      if (_pupper && _pcost) {
1237
        for (int i = 0; i != _arc_num; ++i) {
1238
          Arc e = _arc_ref[i];
1239
          _source[i] = _node_id[_graph.source(e)];
1240
          _target[i] = _node_id[_graph.target(e)];
1241
          _cap[i] = (*_pupper)[e];
1242
          _cost[i] = (*_pcost)[e];
1243
          _flow[i] = 0;
1244
          _state[i] = STATE_LOWER;
1245
        }
1246
      } else {
1247
        for (int i = 0; i != _arc_num; ++i) {
1248
          Arc e = _arc_ref[i];
1249
          _source[i] = _node_id[_graph.source(e)];
1250
          _target[i] = _node_id[_graph.target(e)];
1251
          _flow[i] = 0;
1252
          _state[i] = STATE_LOWER;
1253
        }
1254
        if (_pupper) {
1255
          for (int i = 0; i != _arc_num; ++i)
1256
            _cap[i] = (*_pupper)[_arc_ref[i]];
1257
        } else {
1258
          for (int i = 0; i != _arc_num; ++i)
1259
            _cap[i] = inf_cap;
1260
        }
1261
        if (_pcost) {
1262
          for (int i = 0; i != _arc_num; ++i)
1263
            _cost[i] = (*_pcost)[_arc_ref[i]];
1264
        } else {
1265
          for (int i = 0; i != _arc_num; ++i)
1266
            _cost[i] = 1;
1267
        }
1268
      }
1269
      
1270
      // Remove non-zero lower bounds
1271
      if (_plower) {
1272
        for (int i = 0; i != _arc_num; ++i) {
1273
          Flow c = (*_plower)[_arc_ref[i]];
1274
          if (c != 0) {
1275
            _cap[i] -= c;
1276
            _supply[_source[i]] -= c;
1277
            _supply[_target[i]] += c;
1278
          }
1279
        }
1280
      }
1281

	
1282
      // Add artificial arcs and initialize the spanning tree data structure
1283
      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1284
        _thread[u] = u + 1;
1285
        _rev_thread[u + 1] = u;
1286
        _succ_num[u] = 1;
1287
        _last_succ[u] = u;
1288
        _parent[u] = _root;
1289
        _pred[u] = e;
1290
        _cost[e] = art_cost;
1291
        _cap[e] = inf_cap;
1292
        _state[e] = STATE_TREE;
1293
        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
1294
          _flow[e] = _supply[u];
1295
          _forward[u] = true;
1296
          _pi[u] = -art_cost + _pi[_root];
1297
        } else {
1298
          _flow[e] = -_supply[u];
1299
          _forward[u] = false;
1300
          _pi[u] = art_cost + _pi[_root];
1301
        }
1302
      }
1303

	
1304
      return true;
1305
    }
1306

	
1307
    // Find the join node
1308
    void findJoinNode() {
1309
      int u = _source[in_arc];
1310
      int v = _target[in_arc];
1311
      while (u != v) {
1312
        if (_succ_num[u] < _succ_num[v]) {
1313
          u = _parent[u];
1314
        } else {
1315
          v = _parent[v];
1316
        }
1317
      }
1318
      join = u;
1319
    }
1320

	
1321
    // Find the leaving arc of the cycle and returns true if the
1322
    // leaving arc is not the same as the entering arc
1323
    bool findLeavingArc() {
1324
      // Initialize first and second nodes according to the direction
1325
      // of the cycle
1326
      if (_state[in_arc] == STATE_LOWER) {
1327
        first  = _source[in_arc];
1328
        second = _target[in_arc];
1329
      } else {
1330
        first  = _target[in_arc];
1331
        second = _source[in_arc];
1332
      }
1333
      delta = _cap[in_arc];
1334
      int result = 0;
1335
      Flow d;
1336
      int e;
1337

	
1338
      // Search the cycle along the path form the first node to the root
1339
      for (int u = first; u != join; u = _parent[u]) {
1340
        e = _pred[u];
1341
        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1342
        if (d < delta) {
1343
          delta = d;
1344
          u_out = u;
1345
          result = 1;
1346
        }
1347
      }
1348
      // Search the cycle along the path form the second node to the root
1349
      for (int u = second; u != join; u = _parent[u]) {
1350
        e = _pred[u];
1351
        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1352
        if (d <= delta) {
1353
          delta = d;
1354
          u_out = u;
1355
          result = 2;
1356
        }
1357
      }
1358

	
1359
      if (result == 1) {
1360
        u_in = first;
1361
        v_in = second;
1362
      } else {
1363
        u_in = second;
1364
        v_in = first;
1365
      }
1366
      return result != 0;
1367
    }
1368

	
1369
    // Change _flow and _state vectors
1370
    void changeFlow(bool change) {
1371
      // Augment along the cycle
1372
      if (delta > 0) {
1373
        Flow val = _state[in_arc] * delta;
1374
        _flow[in_arc] += val;
1375
        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1376
          _flow[_pred[u]] += _forward[u] ? -val : val;
1377
        }
1378
        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1379
          _flow[_pred[u]] += _forward[u] ? val : -val;
1380
        }
1381
      }
1382
      // Update the state of the entering and leaving arcs
1383
      if (change) {
1384
        _state[in_arc] = STATE_TREE;
1385
        _state[_pred[u_out]] =
1386
          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1387
      } else {
1388
        _state[in_arc] = -_state[in_arc];
1389
      }
1390
    }
1391

	
1392
    // Update the tree structure
1393
    void updateTreeStructure() {
1394
      int u, w;
1395
      int old_rev_thread = _rev_thread[u_out];
1396
      int old_succ_num = _succ_num[u_out];
1397
      int old_last_succ = _last_succ[u_out];
1398
      v_out = _parent[u_out];
1399

	
1400
      u = _last_succ[u_in];  // the last successor of u_in
1401
      right = _thread[u];    // the node after it
1402

	
1403
      // Handle the case when old_rev_thread equals to v_in
1404
      // (it also means that join and v_out coincide)
1405
      if (old_rev_thread == v_in) {
1406
        last = _thread[_last_succ[u_out]];
1407
      } else {
1408
        last = _thread[v_in];
1409
      }
1410

	
1411
      // Update _thread and _parent along the stem nodes (i.e. the nodes
1412
      // between u_in and u_out, whose parent have to be changed)
1413
      _thread[v_in] = stem = u_in;
1414
      _dirty_revs.clear();
1415
      _dirty_revs.push_back(v_in);
1416
      par_stem = v_in;
1417
      while (stem != u_out) {
1418
        // Insert the next stem node into the thread list
1419
        new_stem = _parent[stem];
1420
        _thread[u] = new_stem;
1421
        _dirty_revs.push_back(u);
1422

	
1423
        // Remove the subtree of stem from the thread list
1424
        w = _rev_thread[stem];
1425
        _thread[w] = right;
1426
        _rev_thread[right] = w;
1427

	
1428
        // Change the parent node and shift stem nodes
1429
        _parent[stem] = par_stem;
1430
        par_stem = stem;
1431
        stem = new_stem;
1432

	
1433
        // Update u and right
1434
        u = _last_succ[stem] == _last_succ[par_stem] ?
1435
          _rev_thread[par_stem] : _last_succ[stem];
1436
        right = _thread[u];
1437
      }
1438
      _parent[u_out] = par_stem;
1439
      _thread[u] = last;
1440
      _rev_thread[last] = u;
1441
      _last_succ[u_out] = u;
1442

	
1443
      // Remove the subtree of u_out from the thread list except for
1444
      // the case when old_rev_thread equals to v_in
1445
      // (it also means that join and v_out coincide)
1446
      if (old_rev_thread != v_in) {
1447
        _thread[old_rev_thread] = right;
1448
        _rev_thread[right] = old_rev_thread;
1449
      }
1450

	
1451
      // Update _rev_thread using the new _thread values
1452
      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1453
        u = _dirty_revs[i];
1454
        _rev_thread[_thread[u]] = u;
1455
      }
1456

	
1457
      // Update _pred, _forward, _last_succ and _succ_num for the
1458
      // stem nodes from u_out to u_in
1459
      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1460
      u = u_out;
1461
      while (u != u_in) {
1462
        w = _parent[u];
1463
        _pred[u] = _pred[w];
1464
        _forward[u] = !_forward[w];
1465
        tmp_sc += _succ_num[u] - _succ_num[w];
1466
        _succ_num[u] = tmp_sc;
1467
        _last_succ[w] = tmp_ls;
1468
        u = w;
1469
      }
1470
      _pred[u_in] = in_arc;
1471
      _forward[u_in] = (u_in == _source[in_arc]);
1472
      _succ_num[u_in] = old_succ_num;
1473

	
1474
      // Set limits for updating _last_succ form v_in and v_out
1475
      // towards the root
1476
      int up_limit_in = -1;
1477
      int up_limit_out = -1;
1478
      if (_last_succ[join] == v_in) {
1479
        up_limit_out = join;
1480
      } else {
1481
        up_limit_in = join;
1482
      }
1483

	
1484
      // Update _last_succ from v_in towards the root
1485
      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1486
           u = _parent[u]) {
1487
        _last_succ[u] = _last_succ[u_out];
1488
      }
1489
      // Update _last_succ from v_out towards the root
1490
      if (join != old_rev_thread && v_in != old_rev_thread) {
1491
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1492
             u = _parent[u]) {
1493
          _last_succ[u] = old_rev_thread;
1494
        }
1495
      } else {
1496
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1497
             u = _parent[u]) {
1498
          _last_succ[u] = _last_succ[u_out];
1499
        }
1500
      }
1501

	
1502
      // Update _succ_num from v_in to join
1503
      for (u = v_in; u != join; u = _parent[u]) {
1504
        _succ_num[u] += old_succ_num;
1505
      }
1506
      // Update _succ_num from v_out to join
1507
      for (u = v_out; u != join; u = _parent[u]) {
1508
        _succ_num[u] -= old_succ_num;
1509
      }
1510
    }
1511

	
1512
    // Update potentials
1513
    void updatePotential() {
1514
      Cost sigma = _forward[u_in] ?
1515
        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1516
        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1517
      // Update potentials in the subtree, which has been moved
1518
      int end = _thread[_last_succ[u_in]];
1519
      for (int u = u_in; u != end; u = _thread[u]) {
1520
        _pi[u] += sigma;
1521
      }
1522
    }
1523

	
1524
    // Execute the algorithm
1525
    bool start(PivotRule pivot_rule) {
1526
      // Select the pivot rule implementation
1527
      switch (pivot_rule) {
1528
        case FIRST_ELIGIBLE:
1529
          return start<FirstEligiblePivotRule>();
1530
        case BEST_ELIGIBLE:
1531
          return start<BestEligiblePivotRule>();
1532
        case BLOCK_SEARCH:
1533
          return start<BlockSearchPivotRule>();
1534
        case CANDIDATE_LIST:
1535
          return start<CandidateListPivotRule>();
1536
        case ALTERING_LIST:
1537
          return start<AlteringListPivotRule>();
1538
      }
1539
      return false;
1540
    }
1541

	
1542
    template <typename PivotRuleImpl>
1543
    bool start() {
1544
      PivotRuleImpl pivot(*this);
1545

	
1546
      // Execute the Network Simplex algorithm
1547
      while (pivot.findEnteringArc()) {
1548
        findJoinNode();
1549
        bool change = findLeavingArc();
1550
        changeFlow(change);
1551
        if (change) {
1552
          updateTreeStructure();
1553
          updatePotential();
1554
        }
1555
      }
1556

	
1557
      // Copy flow values to _flow_map
1558
      if (_plower) {
1559
        for (int i = 0; i != _arc_num; ++i) {
1560
          Arc e = _arc_ref[i];
1561
          _flow_map->set(e, (*_plower)[e] + _flow[i]);
1562
        }
1563
      } else {
1564
        for (int i = 0; i != _arc_num; ++i) {
1565
          _flow_map->set(_arc_ref[i], _flow[i]);
1566
        }
1567
      }
1568
      // Copy potential values to _potential_map
1569
      for (NodeIt n(_graph); n != INVALID; ++n) {
1570
        _potential_map->set(n, _pi[_node_id[n]]);
1571
      }
1572

	
1573
      return true;
1574
    }
1575

	
1576
  }; //class NetworkSimplex
1577

	
1578
  ///@}
1579

	
1580
} //namespace lemon
1581

	
1582
#endif //LEMON_NETWORK_SIMPLEX_H
Ignore white space 6 line context
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

	
22
#include <lemon/list_graph.h>
23
#include <lemon/lgf_reader.h>
24

	
25
#include <lemon/network_simplex.h>
26

	
27
#include <lemon/concepts/digraph.h>
28
#include <lemon/concept_check.h>
29

	
30
#include "test_tools.h"
31

	
32
using namespace lemon;
33

	
34
char test_lgf[] =
35
  "@nodes\n"
36
  "label  sup1 sup2 sup3 sup4 sup5\n"
37
  "    1    20   27    0   20   30\n"
38
  "    2    -4    0    0   -8   -3\n"
39
  "    3     0    0    0    0    0\n"
40
  "    4     0    0    0    0    0\n"
41
  "    5     9    0    0    6   11\n"
42
  "    6    -6    0    0   -5   -6\n"
43
  "    7     0    0    0    0    0\n"
44
  "    8     0    0    0    0    3\n"
45
  "    9     3    0    0    0    0\n"
46
  "   10    -2    0    0   -7   -2\n"
47
  "   11     0    0    0  -10    0\n"
48
  "   12   -20  -27    0  -30  -20\n"
49
  "\n"
50
  "@arcs\n"
51
  "       cost  cap low1 low2\n"
52
  " 1  2    70   11    0    8\n"
53
  " 1  3   150    3    0    1\n"
54
  " 1  4    80   15    0    2\n"
55
  " 2  8    80   12    0    0\n"
56
  " 3  5   140    5    0    3\n"
57
  " 4  6    60   10    0    1\n"
58
  " 4  7    80    2    0    0\n"
59
  " 4  8   110    3    0    0\n"
60
  " 5  7    60   14    0    0\n"
61
  " 5 11   120   12    0    0\n"
62
  " 6  3     0    3    0    0\n"
63
  " 6  9   140    4    0    0\n"
64
  " 6 10    90    8    0    0\n"
65
  " 7  1    30    5    0    0\n"
66
  " 8 12    60   16    0    4\n"
67
  " 9 12    50    6    0    0\n"
68
  "10 12    70   13    0    5\n"
69
  "10  2   100    7    0    0\n"
70
  "10  7    60   10    0    0\n"
71
  "11 10    20   14    0    6\n"
72
  "12 11    30   10    0    0\n"
73
  "\n"
74
  "@attributes\n"
75
  "source 1\n"
76
  "target 12\n";
77

	
78

	
79
enum ProblemType {
80
  EQ,
81
  GEQ,
82
  LEQ
83
};
84

	
85
// Check the interface of an MCF algorithm
86
template <typename GR, typename Flow, typename Cost>
87
class McfClassConcept
88
{
89
public:
90

	
91
  template <typename MCF>
92
  struct Constraints {
93
    void constraints() {
94
      checkConcept<concepts::Digraph, GR>();
95

	
96
      MCF mcf(g);
97

	
98
      b = mcf.reset()
99
             .lowerMap(lower)
100
             .upperMap(upper)
101
             .capacityMap(upper)
102
             .boundMaps(lower, upper)
103
             .costMap(cost)
104
             .supplyMap(sup)
105
             .stSupply(n, n, k)
106
             .flowMap(flow)
107
             .potentialMap(pot)
108
             .run();
109
      
110
      const MCF& const_mcf = mcf;
111

	
112
      const typename MCF::FlowMap &fm = const_mcf.flowMap();
113
      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
114

	
115
      v = const_mcf.totalCost();
116
      double x = const_mcf.template totalCost<double>();
117
      v = const_mcf.flow(a);
118
      v = const_mcf.potential(n);
119

	
120
      ignore_unused_variable_warning(fm);
121
      ignore_unused_variable_warning(pm);
122
      ignore_unused_variable_warning(x);
123
    }
124

	
125
    typedef typename GR::Node Node;
126
    typedef typename GR::Arc Arc;
127
    typedef concepts::ReadMap<Node, Flow> NM;
128
    typedef concepts::ReadMap<Arc, Flow> FAM;
129
    typedef concepts::ReadMap<Arc, Cost> CAM;
130

	
131
    const GR &g;
132
    const FAM &lower;
133
    const FAM &upper;
134
    const CAM &cost;
135
    const NM &sup;
136
    const Node &n;
137
    const Arc &a;
138
    const Flow &k;
139
    Flow v;
140
    bool b;
141

	
142
    typename MCF::FlowMap &flow;
143
    typename MCF::PotentialMap &pot;
144
  };
145

	
146
};
147

	
148

	
149
// Check the feasibility of the given flow (primal soluiton)
150
template < typename GR, typename LM, typename UM,
151
           typename SM, typename FM >
152
bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
153
                const SM& supply, const FM& flow,
154
                ProblemType type = EQ )
155
{
156
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
157

	
158
  for (ArcIt e(gr); e != INVALID; ++e) {
159
    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
160
  }
161

	
162
  for (NodeIt n(gr); n != INVALID; ++n) {
163
    typename SM::Value sum = 0;
164
    for (OutArcIt e(gr, n); e != INVALID; ++e)
165
      sum += flow[e];
166
    for (InArcIt e(gr, n); e != INVALID; ++e)
167
      sum -= flow[e];
168
    bool b = (type ==  EQ && sum == supply[n]) ||
169
             (type == GEQ && sum >= supply[n]) ||
170
             (type == LEQ && sum <= supply[n]);
171
    if (!b) return false;
172
  }
173

	
174
  return true;
175
}
176

	
177
// Check the feasibility of the given potentials (dual soluiton)
178
// using the "Complementary Slackness" optimality condition
179
template < typename GR, typename LM, typename UM,
180
           typename CM, typename SM, typename FM, typename PM >
181
bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
182
                     const CM& cost, const SM& supply, const FM& flow, 
183
                     const PM& pi )
184
{
185
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
186

	
187
  bool opt = true;
188
  for (ArcIt e(gr); opt && e != INVALID; ++e) {
189
    typename CM::Value red_cost =
190
      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
191
    opt = red_cost == 0 ||
192
          (red_cost > 0 && flow[e] == lower[e]) ||
193
          (red_cost < 0 && flow[e] == upper[e]);
194
  }
195
  
196
  for (NodeIt n(gr); opt && n != INVALID; ++n) {
197
    typename SM::Value sum = 0;
198
    for (OutArcIt e(gr, n); e != INVALID; ++e)
199
      sum += flow[e];
200
    for (InArcIt e(gr, n); e != INVALID; ++e)
201
      sum -= flow[e];
202
    opt = (sum == supply[n]) || (pi[n] == 0);
203
  }
204
  
205
  return opt;
206
}
207

	
208
// Run a minimum cost flow algorithm and check the results
209
template < typename MCF, typename GR,
210
           typename LM, typename UM,
211
           typename CM, typename SM >
212
void checkMcf( const MCF& mcf, bool mcf_result,
213
               const GR& gr, const LM& lower, const UM& upper,
214
               const CM& cost, const SM& supply,
215
               bool result, typename CM::Value total,
216
               const std::string &test_id = "",
217
               ProblemType type = EQ )
218
{
219
  check(mcf_result == result, "Wrong result " + test_id);
220
  if (result) {
221
    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
222
          "The flow is not feasible " + test_id);
223
    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
224
    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
225
                         mcf.potentialMap()),
226
          "Wrong potentials " + test_id);
227
  }
228
}
229

	
230
int main()
231
{
232
  // Check the interfaces
233
  {
234
    typedef int Flow;
235
    typedef int Cost;
236
    // TODO: This typedef should be enabled if the standard maps are
237
    // reference maps in the graph concepts (See #190).
238
/**/
239
    //typedef concepts::Digraph GR;
240
    typedef ListDigraph GR;
241
/**/
242
    checkConcept< McfClassConcept<GR, Flow, Cost>,
243
                  NetworkSimplex<GR, Flow, Cost> >();
244
  }
245

	
246
  // Run various MCF tests
247
  typedef ListDigraph Digraph;
248
  DIGRAPH_TYPEDEFS(ListDigraph);
249

	
250
  // Read the test digraph
251
  Digraph gr;
252
  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
253
  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
254
  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
255
  Node v, w;
256

	
257
  std::istringstream input(test_lgf);
258
  DigraphReader<Digraph>(gr, input)
259
    .arcMap("cost", c)
260
    .arcMap("cap", u)
261
    .arcMap("low1", l1)
262
    .arcMap("low2", l2)
263
    .nodeMap("sup1", s1)
264
    .nodeMap("sup2", s2)
265
    .nodeMap("sup3", s3)
266
    .nodeMap("sup4", s4)
267
    .nodeMap("sup5", s5)
268
    .node("source", v)
269
    .node("target", w)
270
    .run();
271

	
272
  // A. Test NetworkSimplex with the default pivot rule
273
  {
274
    NetworkSimplex<Digraph> mcf(gr);
275

	
276
    // Check the equality form
277
    mcf.upperMap(u).costMap(c);
278
    checkMcf(mcf, mcf.supplyMap(s1).run(),
279
             gr, l1, u, c, s1, true,  5240, "#A1");
280
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
281
             gr, l1, u, c, s2, true,  7620, "#A2");
282
    mcf.lowerMap(l2);
283
    checkMcf(mcf, mcf.supplyMap(s1).run(),
284
             gr, l2, u, c, s1, true,  5970, "#A3");
285
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
286
             gr, l2, u, c, s2, true,  8010, "#A4");
287
    mcf.reset();
288
    checkMcf(mcf, mcf.supplyMap(s1).run(),
289
             gr, l1, cu, cc, s1, true,  74, "#A5");
290
    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
291
             gr, l2, cu, cc, s2, true,  94, "#A6");
292
    mcf.reset();
293
    checkMcf(mcf, mcf.run(),
294
             gr, l1, cu, cc, s3, true,   0, "#A7");
295
    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
296
             gr, l2, u, cc, s3, false,   0, "#A8");
297

	
298
    // Check the GEQ form
299
    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
300
    checkMcf(mcf, mcf.run(),
301
             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
302
    mcf.problemType(mcf.GEQ);
303
    checkMcf(mcf, mcf.lowerMap(l2).run(),
304
             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
305
    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
306
    checkMcf(mcf, mcf.run(),
307
             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
308

	
309
    // Check the LEQ form
310
    mcf.reset().problemType(mcf.LEQ);
311
    mcf.upperMap(u).costMap(c).supplyMap(s5);
312
    checkMcf(mcf, mcf.run(),
313
             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
314
    checkMcf(mcf, mcf.lowerMap(l2).run(),
315
             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
316
    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
317
    checkMcf(mcf, mcf.run(),
318
             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
319
  }
320

	
321
  // B. Test NetworkSimplex with each pivot rule
322
  {
323
    NetworkSimplex<Digraph> mcf(gr);
324
    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
325

	
326
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
327
             gr, l2, u, c, s1, true,  5970, "#B1");
328
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
329
             gr, l2, u, c, s1, true,  5970, "#B2");
330
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
331
             gr, l2, u, c, s1, true,  5970, "#B3");
332
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
333
             gr, l2, u, c, s1, true,  5970, "#B4");
334
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
335
             gr, l2, u, c, s1, true,  5970, "#B5");
336
  }
337

	
338
  return 0;
339
}
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
namespace lemon {
20 20

	
21 21
/**
22 22
@defgroup datas Data Structures
23 23
This group contains the several data structures implemented in LEMON.
24 24
*/
25 25

	
26 26
/**
27 27
@defgroup graphs Graph Structures
28 28
@ingroup datas
29 29
\brief Graph structures implemented in LEMON.
30 30

	
31 31
The implementation of combinatorial algorithms heavily relies on
32 32
efficient graph implementations. LEMON offers data structures which are
33 33
planned to be easily used in an experimental phase of implementation studies,
34 34
and thereafter the program code can be made efficient by small modifications.
35 35

	
36 36
The most efficient implementation of diverse applications require the
37 37
usage of different physical graph implementations. These differences
38 38
appear in the size of graph we require to handle, memory or time usage
39 39
limitations or in the set of operations through which the graph can be
40 40
accessed.  LEMON provides several physical graph structures to meet
41 41
the diverging requirements of the possible users.  In order to save on
42 42
running time or on memory usage, some structures may fail to provide
43 43
some graph features like arc/edge or node deletion.
44 44

	
45 45
Alteration of standard containers need a very limited number of
46 46
operations, these together satisfy the everyday requirements.
47 47
In the case of graph structures, different operations are needed which do
48 48
not alter the physical graph, but gives another view. If some nodes or
49 49
arcs have to be hidden or the reverse oriented graph have to be used, then
50 50
this is the case. It also may happen that in a flow implementation
51 51
the residual graph can be accessed by another algorithm, or a node-set
52 52
is to be shrunk for another algorithm.
53 53
LEMON also provides a variety of graphs for these requirements called
54 54
\ref graph_adaptors "graph adaptors". Adaptors cannot be used alone but only
55 55
in conjunction with other graph representations.
56 56

	
57 57
You are free to use the graph structure that fit your requirements
58 58
the best, most graph algorithms and auxiliary data structures can be used
59 59
with any graph structure.
60 60

	
61 61
<b>See also:</b> \ref graph_concepts "Graph Structure Concepts".
62 62
*/
63 63

	
64 64
/**
65 65
@defgroup graph_adaptors Adaptor Classes for Graphs
66 66
@ingroup graphs
67 67
\brief Adaptor classes for digraphs and graphs
68 68

	
69 69
This group contains several useful adaptor classes for digraphs and graphs.
70 70

	
71 71
The main parts of LEMON are the different graph structures, generic
72 72
graph algorithms, graph concepts, which couple them, and graph
73 73
adaptors. While the previous notions are more or less clear, the
74 74
latter one needs further explanation. Graph adaptors are graph classes
75 75
which serve for considering graph structures in different ways.
76 76

	
77 77
A short example makes this much clearer.  Suppose that we have an
78 78
instance \c g of a directed graph type, say ListDigraph and an algorithm
79 79
\code
80 80
template <typename Digraph>
81 81
int algorithm(const Digraph&);
82 82
\endcode
83 83
is needed to run on the reverse oriented graph.  It may be expensive
84 84
(in time or in memory usage) to copy \c g with the reversed
85 85
arcs.  In this case, an adaptor class is used, which (according
86 86
to LEMON \ref concepts::Digraph "digraph concepts") works as a digraph.
87 87
The adaptor uses the original digraph structure and digraph operations when
88 88
methods of the reversed oriented graph are called.  This means that the adaptor
89 89
have minor memory usage, and do not perform sophisticated algorithmic
90 90
actions.  The purpose of it is to give a tool for the cases when a
91 91
graph have to be used in a specific alteration.  If this alteration is
92 92
obtained by a usual construction like filtering the node or the arc set or
93 93
considering a new orientation, then an adaptor is worthwhile to use.
94 94
To come back to the reverse oriented graph, in this situation
95 95
\code
96 96
template<typename Digraph> class ReverseDigraph;
97 97
\endcode
98 98
template class can be used. The code looks as follows
99 99
\code
100 100
ListDigraph g;
101 101
ReverseDigraph<ListDigraph> rg(g);
102 102
int result = algorithm(rg);
103 103
\endcode
104 104
During running the algorithm, the original digraph \c g is untouched.
105 105
This techniques give rise to an elegant code, and based on stable
106 106
graph adaptors, complex algorithms can be implemented easily.
107 107

	
108 108
In flow, circulation and matching problems, the residual
109 109
graph is of particular importance. Combining an adaptor implementing
110 110
this with shortest path algorithms or minimum mean cycle algorithms,
111 111
a range of weighted and cardinality optimization algorithms can be
112 112
obtained. For other examples, the interested user is referred to the
113 113
detailed documentation of particular adaptors.
114 114

	
115 115
The behavior of graph adaptors can be very different. Some of them keep
116 116
capabilities of the original graph while in other cases this would be
117 117
meaningless. This means that the concepts that they meet depend
118 118
on the graph adaptor, and the wrapped graph.
119 119
For example, if an arc of a reversed digraph is deleted, this is carried
120 120
out by deleting the corresponding arc of the original digraph, thus the
121 121
adaptor modifies the original digraph.
122 122
However in case of a residual digraph, this operation has no sense.
123 123

	
124 124
Let us stand one more example here to simplify your work.
125 125
ReverseDigraph has constructor
126 126
\code
127 127
ReverseDigraph(Digraph& digraph);
128 128
\endcode
129 129
This means that in a situation, when a <tt>const %ListDigraph&</tt>
130 130
reference to a graph is given, then it have to be instantiated with
131 131
<tt>Digraph=const %ListDigraph</tt>.
132 132
\code
133 133
int algorithm1(const ListDigraph& g) {
134 134
  ReverseDigraph<const ListDigraph> rg(g);
135 135
  return algorithm2(rg);
136 136
}
137 137
\endcode
138 138
*/
139 139

	
140 140
/**
141 141
@defgroup semi_adaptors Semi-Adaptor Classes for Graphs
142 142
@ingroup graphs
143 143
\brief Graph types between real graphs and graph adaptors.
144 144

	
145 145
This group contains some graph types between real graphs and graph adaptors.
146 146
These classes wrap graphs to give new functionality as the adaptors do it.
147 147
On the other hand they are not light-weight structures as the adaptors.
148 148
*/
149 149

	
150 150
/**
151 151
@defgroup maps Maps
152 152
@ingroup datas
153 153
\brief Map structures implemented in LEMON.
154 154

	
155 155
This group contains the map structures implemented in LEMON.
156 156

	
157 157
LEMON provides several special purpose maps and map adaptors that e.g. combine
158 158
new maps from existing ones.
159 159

	
160 160
<b>See also:</b> \ref map_concepts "Map Concepts".
161 161
*/
162 162

	
163 163
/**
164 164
@defgroup graph_maps Graph Maps
165 165
@ingroup maps
166 166
\brief Special graph-related maps.
167 167

	
168 168
This group contains maps that are specifically designed to assign
169 169
values to the nodes and arcs/edges of graphs.
170 170

	
171 171
If you are looking for the standard graph maps (\c NodeMap, \c ArcMap,
172 172
\c EdgeMap), see the \ref graph_concepts "Graph Structure Concepts".
173 173
*/
174 174

	
175 175
/**
176 176
\defgroup map_adaptors Map Adaptors
177 177
\ingroup maps
178 178
\brief Tools to create new maps from existing ones
179 179

	
180 180
This group contains map adaptors that are used to create "implicit"
181 181
maps from other maps.
182 182

	
183 183
Most of them are \ref concepts::ReadMap "read-only maps".
184 184
They can make arithmetic and logical operations between one or two maps
185 185
(negation, shifting, addition, multiplication, logical 'and', 'or',
186 186
'not' etc.) or e.g. convert a map to another one of different Value type.
187 187

	
188 188
The typical usage of this classes is passing implicit maps to
189 189
algorithms.  If a function type algorithm is called then the function
190 190
type map adaptors can be used comfortable. For example let's see the
191 191
usage of map adaptors with the \c graphToEps() function.
192 192
\code
193 193
  Color nodeColor(int deg) {
194 194
    if (deg >= 2) {
195 195
      return Color(0.5, 0.0, 0.5);
196 196
    } else if (deg == 1) {
197 197
      return Color(1.0, 0.5, 1.0);
198 198
    } else {
199 199
      return Color(0.0, 0.0, 0.0);
200 200
    }
201 201
  }
202 202

	
203 203
  Digraph::NodeMap<int> degree_map(graph);
204 204

	
205 205
  graphToEps(graph, "graph.eps")
206 206
    .coords(coords).scaleToA4().undirected()
207 207
    .nodeColors(composeMap(functorToMap(nodeColor), degree_map))
208 208
    .run();
209 209
\endcode
210 210
The \c functorToMap() function makes an \c int to \c Color map from the
211 211
\c nodeColor() function. The \c composeMap() compose the \c degree_map
212 212
and the previously created map. The composed map is a proper function to
213 213
get the color of each node.
214 214

	
215 215
The usage with class type algorithms is little bit harder. In this
216 216
case the function type map adaptors can not be used, because the
217 217
function map adaptors give back temporary objects.
218 218
\code
219 219
  Digraph graph;
220 220

	
221 221
  typedef Digraph::ArcMap<double> DoubleArcMap;
222 222
  DoubleArcMap length(graph);
223 223
  DoubleArcMap speed(graph);
224 224

	
225 225
  typedef DivMap<DoubleArcMap, DoubleArcMap> TimeMap;
226 226
  TimeMap time(length, speed);
227 227

	
228 228
  Dijkstra<Digraph, TimeMap> dijkstra(graph, time);
229 229
  dijkstra.run(source, target);
230 230
\endcode
231 231
We have a length map and a maximum speed map on the arcs of a digraph.
232 232
The minimum time to pass the arc can be calculated as the division of
233 233
the two maps which can be done implicitly with the \c DivMap template
234 234
class. We use the implicit minimum time map as the length map of the
235 235
\c Dijkstra algorithm.
236 236
*/
237 237

	
238 238
/**
239 239
@defgroup matrices Matrices
240 240
@ingroup datas
241 241
\brief Two dimensional data storages implemented in LEMON.
242 242

	
243 243
This group contains two dimensional data storages implemented in LEMON.
244 244
*/
245 245

	
246 246
/**
247 247
@defgroup paths Path Structures
248 248
@ingroup datas
249 249
\brief %Path structures implemented in LEMON.
250 250

	
251 251
This group contains the path structures implemented in LEMON.
252 252

	
253 253
LEMON provides flexible data structures to work with paths.
254 254
All of them have similar interfaces and they can be copied easily with
255 255
assignment operators and copy constructors. This makes it easy and
256 256
efficient to have e.g. the Dijkstra algorithm to store its result in
257 257
any kind of path structure.
258 258

	
259 259
\sa lemon::concepts::Path
260 260
*/
261 261

	
262 262
/**
263 263
@defgroup auxdat Auxiliary Data Structures
264 264
@ingroup datas
265 265
\brief Auxiliary data structures implemented in LEMON.
266 266

	
267 267
This group contains some data structures implemented in LEMON in
268 268
order to make it easier to implement combinatorial algorithms.
269 269
*/
270 270

	
271 271
/**
272 272
@defgroup algs Algorithms
273 273
\brief This group contains the several algorithms
274 274
implemented in LEMON.
275 275

	
276 276
This group contains the several algorithms
277 277
implemented in LEMON.
278 278
*/
279 279

	
280 280
/**
281 281
@defgroup search Graph Search
282 282
@ingroup algs
283 283
\brief Common graph search algorithms.
284 284

	
285 285
This group contains the common graph search algorithms, namely
286 286
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS).
287 287
*/
288 288

	
289 289
/**
290 290
@defgroup shortest_path Shortest Path Algorithms
291 291
@ingroup algs
292 292
\brief Algorithms for finding shortest paths.
293 293

	
294 294
This group contains the algorithms for finding shortest paths in digraphs.
295 295

	
296 296
 - \ref Dijkstra algorithm for finding shortest paths from a source node
297 297
   when all arc lengths are non-negative.
298 298
 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths
299 299
   from a source node when arc lenghts can be either positive or negative,
300 300
   but the digraph should not contain directed cycles with negative total
301 301
   length.
302 302
 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms
303 303
   for solving the \e all-pairs \e shortest \e paths \e problem when arc
304 304
   lenghts can be either positive or negative, but the digraph should
305 305
   not contain directed cycles with negative total length.
306 306
 - \ref Suurballe A successive shortest path algorithm for finding
307 307
   arc-disjoint paths between two nodes having minimum total length.
308 308
*/
309 309

	
310 310
/**
311 311
@defgroup max_flow Maximum Flow Algorithms
312 312
@ingroup algs
313 313
\brief Algorithms for finding maximum flows.
314 314

	
315 315
This group contains the algorithms for finding maximum flows and
316 316
feasible circulations.
317 317

	
318 318
The \e maximum \e flow \e problem is to find a flow of maximum value between
319 319
a single source and a single target. Formally, there is a \f$G=(V,A)\f$
320
digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
320
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
321 321
\f$s, t \in V\f$ source and target nodes.
322
A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
322
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
323 323
following optimization problem.
324 324

	
325
\f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
326
\f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
327
    \qquad \forall v\in V\setminus\{s,t\} \f]
328
\f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
325
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
326
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
327
    \quad \forall u\in V\setminus\{s,t\} \f]
328
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
329 329

	
330 330
LEMON contains several algorithms for solving maximum flow problems:
331 331
- \ref EdmondsKarp Edmonds-Karp algorithm.
332 332
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm.
333 333
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees.
334 334
- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees.
335 335

	
336 336
In most cases the \ref Preflow "Preflow" algorithm provides the
337 337
fastest method for computing a maximum flow. All implementations
338 338
provides functions to also query the minimum cut, which is the dual
339 339
problem of the maximum flow.
340 340
*/
341 341

	
342 342
/**
343 343
@defgroup min_cost_flow Minimum Cost Flow Algorithms
344 344
@ingroup algs
345 345

	
346 346
\brief Algorithms for finding minimum cost flows and circulations.
347 347

	
348 348
This group contains the algorithms for finding minimum cost flows and
349 349
circulations.
350 350

	
351 351
The \e minimum \e cost \e flow \e problem is to find a feasible flow of
352 352
minimum total cost from a set of supply nodes to a set of demand nodes
353
in a network with capacity constraints and arc costs.
353
in a network with capacity constraints (lower and upper bounds)
354
and arc costs.
354 355
Formally, let \f$G=(V,A)\f$ be a digraph,
355 356
\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
356
upper bounds for the flow values on the arcs,
357
upper bounds for the flow values on the arcs, for which
358
\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
357 359
\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
358
on the arcs, and
359
\f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
360
of the nodes.
361
A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
362
the following optimization problem.
360
on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
361
signed supply values of the nodes.
362
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
363
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
364
\f$-sup(u)\f$ demand.
365
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
366
of the following optimization problem.
363 367

	
364
\f[ \min\sum_{a\in A} f(a) cost(a) \f]
365
\f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
366
    supply(v) \qquad \forall v\in V \f]
367
\f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
368
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
369
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
370
    sup(u) \quad \forall u\in V \f]
371
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
368 372

	
369
LEMON contains several algorithms for solving minimum cost flow problems:
370
 - \ref CycleCanceling Cycle-canceling algorithms.
371
 - \ref CapacityScaling Successive shortest path algorithm with optional
373
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
374
zero or negative in order to have a feasible solution (since the sum
375
of the expressions on the left-hand side of the inequalities is zero).
376
It means that the total demand must be greater or equal to the total
377
supply and all the supplies have to be carried out from the supply nodes,
378
but there could be demands that are not satisfied.
379
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
380
constraints have to be satisfied with equality, i.e. all demands
381
have to be satisfied and all supplies have to be used.
382

	
383
If you need the opposite inequalities in the supply/demand constraints
384
(i.e. the total demand is less than the total supply and all the demands
385
have to be satisfied while there could be supplies that are not used),
386
then you could easily transform the problem to the above form by reversing
387
the direction of the arcs and taking the negative of the supply values
388
(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
389
However \ref NetworkSimplex algorithm also supports this form directly
390
for the sake of convenience.
391

	
392
A feasible solution for this problem can be found using \ref Circulation.
393

	
394
Note that the above formulation is actually more general than the usual
395
definition of the minimum cost flow problem, in which strict equalities
396
are required in the supply/demand contraints, i.e.
397

	
398
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
399
    sup(u) \quad \forall u\in V. \f]
400

	
401
However if the sum of the supply values is zero, then these two problems
402
are equivalent. So if you need the equality form, you have to ensure this
403
additional contraint for the algorithms.
404

	
405
The dual solution of the minimum cost flow problem is represented by node 
406
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
407
An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
408
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
409
node potentials the following \e complementary \e slackness optimality
410
conditions hold.
411

	
412
 - For all \f$uv\in A\f$ arcs:
413
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
414
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
415
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
416
 - For all \f$u\in V\f$:
417
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
418
     then \f$\pi(u)=0\f$.
419
 
420
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
421
\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
422
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
423

	
424
All algorithms provide dual solution (node potentials) as well
425
if an optimal flow is found.
426

	
427
LEMON contains several algorithms for solving minimum cost flow problems.
428
 - \ref NetworkSimplex Primal Network Simplex algorithm with various
429
   pivot strategies.
430
 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
431
   cost scaling.
432
 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
372 433
   capacity scaling.
373
 - \ref CostScaling Push-relabel and augment-relabel algorithms based on
374
   cost scaling.
375
 - \ref NetworkSimplex Primal network simplex algorithm with various
376
   pivot strategies.
434
 - \ref CancelAndTighten The Cancel and Tighten algorithm.
435
 - \ref CycleCanceling Cycle-Canceling algorithms.
436

	
437
Most of these implementations support the general inequality form of the
438
minimum cost flow problem, but CancelAndTighten and CycleCanceling
439
only support the equality form due to the primal method they use.
440

	
441
In general NetworkSimplex is the most efficient implementation,
442
but in special cases other algorithms could be faster.
443
For example, if the total supply and/or capacities are rather small,
444
CapacityScaling is usually the fastest algorithm (without effective scaling).
377 445
*/
378 446

	
379 447
/**
380 448
@defgroup min_cut Minimum Cut Algorithms
381 449
@ingroup algs
382 450

	
383 451
\brief Algorithms for finding minimum cut in graphs.
384 452

	
385 453
This group contains the algorithms for finding minimum cut in graphs.
386 454

	
387 455
The \e minimum \e cut \e problem is to find a non-empty and non-complete
388 456
\f$X\f$ subset of the nodes with minimum overall capacity on
389 457
outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a
390 458
\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum
391 459
cut is the \f$X\f$ solution of the next optimization problem:
392 460

	
393 461
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}
394 462
    \sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f]
395 463

	
396 464
LEMON contains several algorithms related to minimum cut problems:
397 465

	
398 466
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
399 467
  in directed graphs.
400 468
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
401 469
  calculating minimum cut in undirected graphs.
402 470
- \ref GomoryHu "Gomory-Hu tree computation" for calculating
403 471
  all-pairs minimum cut in undirected graphs.
404 472

	
405 473
If you want to find minimum cut just between two distinict nodes,
406 474
see the \ref max_flow "maximum flow problem".
407 475
*/
408 476

	
409 477
/**
410 478
@defgroup graph_properties Connectivity and Other Graph Properties
411 479
@ingroup algs
412 480
\brief Algorithms for discovering the graph properties
413 481

	
414 482
This group contains the algorithms for discovering the graph properties
415 483
like connectivity, bipartiteness, euler property, simplicity etc.
416 484

	
417 485
\image html edge_biconnected_components.png
418 486
\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth
419 487
*/
420 488

	
421 489
/**
422 490
@defgroup planar Planarity Embedding and Drawing
423 491
@ingroup algs
424 492
\brief Algorithms for planarity checking, embedding and drawing
425 493

	
426 494
This group contains the algorithms for planarity checking,
427 495
embedding and drawing.
428 496

	
429 497
\image html planar.png
430 498
\image latex planar.eps "Plane graph" width=\textwidth
431 499
*/
432 500

	
433 501
/**
434 502
@defgroup matching Matching Algorithms
435 503
@ingroup algs
436 504
\brief Algorithms for finding matchings in graphs and bipartite graphs.
437 505

	
438 506
This group contains the algorithms for calculating
439 507
matchings in graphs and bipartite graphs. The general matching problem is
440 508
finding a subset of the edges for which each node has at most one incident
441 509
edge.
442 510

	
443 511
There are several different algorithms for calculate matchings in
444 512
graphs.  The matching problems in bipartite graphs are generally
445 513
easier than in general graphs. The goal of the matching optimization
446 514
can be finding maximum cardinality, maximum weight or minimum cost
447 515
matching. The search can be constrained to find perfect or
448 516
maximum cardinality matching.
449 517

	
450 518
The matching algorithms implemented in LEMON:
451 519
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
452 520
  for calculating maximum cardinality matching in bipartite graphs.
453 521
- \ref PrBipartiteMatching Push-relabel algorithm
454 522
  for calculating maximum cardinality matching in bipartite graphs.
455 523
- \ref MaxWeightedBipartiteMatching
456 524
  Successive shortest path algorithm for calculating maximum weighted
457 525
  matching and maximum weighted bipartite matching in bipartite graphs.
458 526
- \ref MinCostMaxBipartiteMatching
459 527
  Successive shortest path algorithm for calculating minimum cost maximum
460 528
  matching in bipartite graphs.
461 529
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
462 530
  maximum cardinality matching in general graphs.
463 531
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
464 532
  maximum weighted matching in general graphs.
465 533
- \ref MaxWeightedPerfectMatching
466 534
  Edmond's blossom shrinking algorithm for calculating maximum weighted
467 535
  perfect matching in general graphs.
468 536

	
469 537
\image html bipartite_matching.png
470 538
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
471 539
*/
472 540

	
473 541
/**
474 542
@defgroup spantree Minimum Spanning Tree Algorithms
475 543
@ingroup algs
476 544
\brief Algorithms for finding a minimum cost spanning tree in a graph.
477 545

	
478 546
This group contains the algorithms for finding a minimum cost spanning
479 547
tree in a graph.
480 548
*/
481 549

	
482 550
/**
483 551
@defgroup auxalg Auxiliary Algorithms
484 552
@ingroup algs
485 553
\brief Auxiliary algorithms implemented in LEMON.
486 554

	
487 555
This group contains some algorithms implemented in LEMON
488 556
in order to make it easier to implement complex algorithms.
489 557
*/
490 558

	
491 559
/**
492 560
@defgroup approx Approximation Algorithms
493 561
@ingroup algs
494 562
\brief Approximation algorithms.
495 563

	
496 564
This group contains the approximation and heuristic algorithms
497 565
implemented in LEMON.
498 566
*/
499 567

	
500 568
/**
501 569
@defgroup gen_opt_group General Optimization Tools
502 570
\brief This group contains some general optimization frameworks
503 571
implemented in LEMON.
504 572

	
505 573
This group contains some general optimization frameworks
506 574
implemented in LEMON.
507 575
*/
508 576

	
509 577
/**
510 578
@defgroup lp_group Lp and Mip Solvers
511 579
@ingroup gen_opt_group
512 580
\brief Lp and Mip solver interfaces for LEMON.
513 581

	
514 582
This group contains Lp and Mip solver interfaces for LEMON. The
515 583
various LP solvers could be used in the same manner with this
516 584
interface.
517 585
*/
518 586

	
519 587
/**
520 588
@defgroup lp_utils Tools for Lp and Mip Solvers
521 589
@ingroup lp_group
522 590
\brief Helper tools to the Lp and Mip solvers.
523 591

	
524 592
This group adds some helper tools to general optimization framework
525 593
implemented in LEMON.
526 594
*/
527 595

	
528 596
/**
529 597
@defgroup metah Metaheuristics
530 598
@ingroup gen_opt_group
531 599
\brief Metaheuristics for LEMON library.
532 600

	
533 601
This group contains some metaheuristic optimization tools.
534 602
*/
535 603

	
536 604
/**
537 605
@defgroup utils Tools and Utilities
538 606
\brief Tools and utilities for programming in LEMON
539 607

	
540 608
Tools and utilities for programming in LEMON.
541 609
*/
542 610

	
543 611
/**
544 612
@defgroup gutils Basic Graph Utilities
545 613
@ingroup utils
546 614
\brief Simple basic graph utilities.
547 615

	
548 616
This group contains some simple basic graph utilities.
549 617
*/
550 618

	
551 619
/**
552 620
@defgroup misc Miscellaneous Tools
553 621
@ingroup utils
554 622
\brief Tools for development, debugging and testing.
555 623

	
556 624
This group contains several useful tools for development,
557 625
debugging and testing.
558 626
*/
559 627

	
560 628
/**
561 629
@defgroup timecount Time Measuring and Counting
562 630
@ingroup misc
563 631
\brief Simple tools for measuring the performance of algorithms.
564 632

	
565 633
This group contains simple tools for measuring the performance
566 634
of algorithms.
567 635
*/
568 636

	
569 637
/**
570 638
@defgroup exceptions Exceptions
571 639
@ingroup utils
572 640
\brief Exceptions defined in LEMON.
573 641

	
574 642
This group contains the exceptions defined in LEMON.
575 643
*/
576 644

	
577 645
/**
578 646
@defgroup io_group Input-Output
579 647
\brief Graph Input-Output methods
580 648

	
581 649
This group contains the tools for importing and exporting graphs
582 650
and graph related data. Now it supports the \ref lgf-format
583 651
"LEMON Graph Format", the \c DIMACS format and the encapsulated
584 652
postscript (EPS) format.
585 653
*/
586 654

	
587 655
/**
588 656
@defgroup lemon_io LEMON Graph Format
589 657
@ingroup io_group
590 658
\brief Reading and writing LEMON Graph Format.
591 659

	
592 660
This group contains methods for reading and writing
593 661
\ref lgf-format "LEMON Graph Format".
594 662
*/
595 663

	
596 664
/**
597 665
@defgroup eps_io Postscript Exporting
598 666
@ingroup io_group
599 667
\brief General \c EPS drawer and graph exporter
600 668

	
601 669
This group contains general \c EPS drawing methods and special
602 670
graph exporting tools.
603 671
*/
604 672

	
605 673
/**
606 674
@defgroup dimacs_group DIMACS format
607 675
@ingroup io_group
608 676
\brief Read and write files in DIMACS format
609 677

	
610 678
Tools to read a digraph from or write it to a file in DIMACS format data.
611 679
*/
612 680

	
613 681
/**
614 682
@defgroup nauty_group NAUTY Format
615 683
@ingroup io_group
616 684
\brief Read \e Nauty format
617 685

	
618 686
Tool to read graphs from \e Nauty format data.
619 687
*/
620 688

	
621 689
/**
622 690
@defgroup concept Concepts
623 691
\brief Skeleton classes and concept checking classes
624 692

	
625 693
This group contains the data/algorithm skeletons and concept checking
626 694
classes implemented in LEMON.
627 695

	
628 696
The purpose of the classes in this group is fourfold.
629 697

	
630 698
- These classes contain the documentations of the %concepts. In order
631 699
  to avoid document multiplications, an implementation of a concept
632 700
  simply refers to the corresponding concept class.
633 701

	
634 702
- These classes declare every functions, <tt>typedef</tt>s etc. an
635 703
  implementation of the %concepts should provide, however completely
636 704
  without implementations and real data structures behind the
637 705
  interface. On the other hand they should provide nothing else. All
638 706
  the algorithms working on a data structure meeting a certain concept
639 707
  should compile with these classes. (Though it will not run properly,
640 708
  of course.) In this way it is easily to check if an algorithm
641 709
  doesn't use any extra feature of a certain implementation.
642 710

	
643 711
- The concept descriptor classes also provide a <em>checker class</em>
644 712
  that makes it possible to check whether a certain implementation of a
645 713
  concept indeed provides all the required features.
646 714

	
647 715
- Finally, They can serve as a skeleton of a new implementation of a concept.
648 716
*/
649 717

	
650 718
/**
651 719
@defgroup graph_concepts Graph Structure Concepts
652 720
@ingroup concept
653 721
\brief Skeleton and concept checking classes for graph structures
654 722

	
655 723
This group contains the skeletons and concept checking classes of LEMON's
656 724
graph structures and helper classes used to implement these.
657 725
*/
658 726

	
659 727
/**
660 728
@defgroup map_concepts Map Concepts
661 729
@ingroup concept
662 730
\brief Skeleton and concept checking classes for maps
663 731

	
664 732
This group contains the skeletons and concept checking classes of maps.
665 733
*/
666 734

	
667 735
/**
668 736
\anchor demoprograms
669 737

	
670 738
@defgroup demos Demo Programs
671 739

	
672 740
Some demo programs are listed here. Their full source codes can be found in
673 741
the \c demo subdirectory of the source tree.
674 742

	
675 743
In order to compile them, use the <tt>make demo</tt> or the
676 744
<tt>make check</tt> commands.
677 745
*/
678 746

	
679 747
/**
680 748
@defgroup tools Standalone Utility Applications
681 749

	
682 750
Some utility applications are listed here.
683 751

	
684 752
The standard compilation procedure (<tt>./configure;make</tt>) will compile
685 753
them, as well.
686 754
*/
687 755

	
688 756
}
Ignore white space 6 line context
1 1
EXTRA_DIST += \
2 2
	lemon/lemon.pc.in \
3 3
	lemon/CMakeLists.txt
4 4

	
5 5
pkgconfig_DATA += lemon/lemon.pc
6 6

	
7 7
lib_LTLIBRARIES += lemon/libemon.la
8 8

	
9 9
lemon_libemon_la_SOURCES = \
10 10
	lemon/arg_parser.cc \
11 11
	lemon/base.cc \
12 12
	lemon/color.cc \
13 13
	lemon/lp_base.cc \
14 14
	lemon/lp_skeleton.cc \
15 15
	lemon/random.cc \
16 16
	lemon/bits/windows.cc
17 17

	
18 18

	
19 19
lemon_libemon_la_CXXFLAGS = \
20 20
	$(AM_CXXFLAGS) \
21 21
	$(GLPK_CFLAGS) \
22 22
	$(CPLEX_CFLAGS) \
23 23
	$(SOPLEX_CXXFLAGS) \
24 24
	$(CLP_CXXFLAGS) \
25 25
	$(CBC_CXXFLAGS)
26 26

	
27 27
lemon_libemon_la_LDFLAGS = \
28 28
	$(GLPK_LIBS) \
29 29
	$(CPLEX_LIBS) \
30 30
	$(SOPLEX_LIBS) \
31 31
	$(CLP_LIBS) \
32 32
	$(CBC_LIBS)
33 33

	
34 34
if HAVE_GLPK
35 35
lemon_libemon_la_SOURCES += lemon/glpk.cc
36 36
endif
37 37

	
38 38
if HAVE_CPLEX
39 39
lemon_libemon_la_SOURCES += lemon/cplex.cc
40 40
endif
41 41

	
42 42
if HAVE_SOPLEX
43 43
lemon_libemon_la_SOURCES += lemon/soplex.cc
44 44
endif
45 45

	
46 46
if HAVE_CLP
47 47
lemon_libemon_la_SOURCES += lemon/clp.cc
48 48
endif
49 49

	
50 50
if HAVE_CBC
51 51
lemon_libemon_la_SOURCES += lemon/cbc.cc
52 52
endif
53 53

	
54 54
lemon_HEADERS += \
55 55
	lemon/adaptors.h \
56 56
	lemon/arg_parser.h \
57 57
	lemon/assert.h \
58 58
	lemon/bfs.h \
59 59
	lemon/bin_heap.h \
60 60
	lemon/circulation.h \
61 61
	lemon/clp.h \
62 62
	lemon/color.h \
63 63
	lemon/concept_check.h \
64 64
	lemon/connectivity.h \
65 65
	lemon/counter.h \
66 66
	lemon/core.h \
67 67
	lemon/cplex.h \
68 68
	lemon/dfs.h \
69 69
	lemon/dijkstra.h \
70 70
	lemon/dim2.h \
71 71
	lemon/dimacs.h \
72 72
	lemon/edge_set.h \
73 73
	lemon/elevator.h \
74 74
	lemon/error.h \
75 75
	lemon/euler.h \
76 76
	lemon/full_graph.h \
77 77
	lemon/glpk.h \
78 78
	lemon/gomory_hu.h \
79 79
	lemon/graph_to_eps.h \
80 80
	lemon/grid_graph.h \
81 81
	lemon/hypercube_graph.h \
82 82
	lemon/kruskal.h \
83 83
	lemon/hao_orlin.h \
84 84
	lemon/lgf_reader.h \
85 85
	lemon/lgf_writer.h \
86 86
	lemon/list_graph.h \
87 87
	lemon/lp.h \
88 88
	lemon/lp_base.h \
89 89
	lemon/lp_skeleton.h \
90 90
	lemon/list_graph.h \
91 91
	lemon/maps.h \
92 92
	lemon/matching.h \
93 93
	lemon/math.h \
94 94
	lemon/min_cost_arborescence.h \
95 95
	lemon/nauty_reader.h \
96
	lemon/network_simplex.h \
96 97
	lemon/path.h \
97 98
	lemon/preflow.h \
98 99
	lemon/radix_sort.h \
99 100
	lemon/random.h \
100 101
	lemon/smart_graph.h \
101 102
	lemon/soplex.h \
102 103
	lemon/suurballe.h \
103 104
	lemon/time_measure.h \
104 105
	lemon/tolerance.h \
105 106
	lemon/unionfind.h \
106 107
	lemon/bits/windows.h
107 108

	
108 109
bits_HEADERS += \
109 110
	lemon/bits/alteration_notifier.h \
110 111
	lemon/bits/array_map.h \
111 112
	lemon/bits/base_extender.h \
112 113
	lemon/bits/bezier.h \
113 114
	lemon/bits/default_map.h \
114 115
	lemon/bits/edge_set_extender.h \
115 116
	lemon/bits/enable_if.h \
116 117
	lemon/bits/graph_adaptor_extender.h \
117 118
	lemon/bits/graph_extender.h \
118 119
	lemon/bits/map_extender.h \
119 120
	lemon/bits/path_dump.h \
120 121
	lemon/bits/solver_bits.h \
121 122
	lemon/bits/traits.h \
122 123
	lemon/bits/variant.h \
123 124
	lemon/bits/vector_map.h
124 125

	
125 126
concept_HEADERS += \
126 127
	lemon/concepts/digraph.h \
127 128
	lemon/concepts/graph.h \
128 129
	lemon/concepts/graph_components.h \
129 130
	lemon/concepts/heap.h \
130 131
	lemon/concepts/maps.h \
131 132
	lemon/concepts/path.h
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_CIRCULATION_H
20 20
#define LEMON_CIRCULATION_H
21 21

	
22 22
#include <lemon/tolerance.h>
23 23
#include <lemon/elevator.h>
24 24

	
25 25
///\ingroup max_flow
26 26
///\file
27 27
///\brief Push-relabel algorithm for finding a feasible circulation.
28 28
///
29 29
namespace lemon {
30 30

	
31 31
  /// \brief Default traits class of Circulation class.
32 32
  ///
33 33
  /// Default traits class of Circulation class.
34
  /// \tparam GR Digraph type.
35
  /// \tparam LM Lower bound capacity map type.
36
  /// \tparam UM Upper bound capacity map type.
37
  /// \tparam DM Delta map type.
34
  ///
35
  /// \tparam GR Type of the digraph the algorithm runs on.
36
  /// \tparam LM The type of the lower bound map.
37
  /// \tparam UM The type of the upper bound (capacity) map.
38
  /// \tparam SM The type of the supply map.
38 39
  template <typename GR, typename LM,
39
            typename UM, typename DM>
40
            typename UM, typename SM>
40 41
  struct CirculationDefaultTraits {
41 42

	
42 43
    /// \brief The type of the digraph the algorithm runs on.
43 44
    typedef GR Digraph;
44 45

	
45
    /// \brief The type of the map that stores the circulation lower
46
    /// bound.
46
    /// \brief The type of the lower bound map.
47 47
    ///
48
    /// The type of the map that stores the circulation lower bound.
49
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
50
    typedef LM LCapMap;
48
    /// The type of the map that stores the lower bounds on the arcs.
49
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
50
    typedef LM LowerMap;
51 51

	
52
    /// \brief The type of the map that stores the circulation upper
53
    /// bound.
52
    /// \brief The type of the upper bound (capacity) map.
54 53
    ///
55
    /// The type of the map that stores the circulation upper bound.
56
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
57
    typedef UM UCapMap;
54
    /// The type of the map that stores the upper bounds (capacities)
55
    /// on the arcs.
56
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
57
    typedef UM UpperMap;
58 58

	
59
    /// \brief The type of the map that stores the lower bound for
60
    /// the supply of the nodes.
59
    /// \brief The type of supply map.
61 60
    ///
62
    /// The type of the map that stores the lower bound for the supply
63
    /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
64
    /// concept.
65
    typedef DM DeltaMap;
61
    /// The type of the map that stores the signed supply values of the 
62
    /// nodes. 
63
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
64
    typedef SM SupplyMap;
66 65

	
67 66
    /// \brief The type of the flow values.
68
    typedef typename DeltaMap::Value Value;
67
    typedef typename SupplyMap::Value Flow;
69 68

	
70 69
    /// \brief The type of the map that stores the flow values.
71 70
    ///
72 71
    /// The type of the map that stores the flow values.
73
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
74
    typedef typename Digraph::template ArcMap<Value> FlowMap;
72
    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
73
    /// concept.
74
    typedef typename Digraph::template ArcMap<Flow> FlowMap;
75 75

	
76 76
    /// \brief Instantiates a FlowMap.
77 77
    ///
78 78
    /// This function instantiates a \ref FlowMap.
79
    /// \param digraph The digraph, to which we would like to define
79
    /// \param digraph The digraph for which we would like to define
80 80
    /// the flow map.
81 81
    static FlowMap* createFlowMap(const Digraph& digraph) {
82 82
      return new FlowMap(digraph);
83 83
    }
84 84

	
85 85
    /// \brief The elevator type used by the algorithm.
86 86
    ///
87 87
    /// The elevator type used by the algorithm.
88 88
    ///
89 89
    /// \sa Elevator
90 90
    /// \sa LinkedElevator
91 91
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
92 92

	
93 93
    /// \brief Instantiates an Elevator.
94 94
    ///
95 95
    /// This function instantiates an \ref Elevator.
96
    /// \param digraph The digraph, to which we would like to define
96
    /// \param digraph The digraph for which we would like to define
97 97
    /// the elevator.
98 98
    /// \param max_level The maximum level of the elevator.
99 99
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
100 100
      return new Elevator(digraph, max_level);
101 101
    }
102 102

	
103 103
    /// \brief The tolerance used by the algorithm
104 104
    ///
105 105
    /// The tolerance used by the algorithm to handle inexact computation.
106
    typedef lemon::Tolerance<Value> Tolerance;
106
    typedef lemon::Tolerance<Flow> Tolerance;
107 107

	
108 108
  };
109 109

	
110 110
  /**
111 111
     \brief Push-relabel algorithm for the network circulation problem.
112 112

	
113 113
     \ingroup max_flow
114
     This class implements a push-relabel algorithm for the network
115
     circulation problem.
114
     This class implements a push-relabel algorithm for the \e network
115
     \e circulation problem.
116 116
     It is to find a feasible circulation when lower and upper bounds
117
     are given for the flow values on the arcs and lower bounds
118
     are given for the supply values of the nodes.
117
     are given for the flow values on the arcs and lower bounds are
118
     given for the difference between the outgoing and incoming flow
119
     at the nodes.
119 120

	
120 121
     The exact formulation of this problem is the following.
121 122
     Let \f$G=(V,A)\f$ be a digraph,
122
     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
123
     \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
124
     \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
125
     \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
126
     \geq delta(v) \quad \forall v\in V, \f]
127
     \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
128
     \note \f$delta(v)\f$ specifies a lower bound for the supply of node
129
     \f$v\f$. It can be either positive or negative, however note that
130
     \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
131
     have a feasible solution.
123
     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$ denote the lower and
124
     upper bounds on the arcs, for which \f$0 \leq lower(uv) \leq upper(uv)\f$
125
     holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
126
     denotes the signed supply values of the nodes.
127
     If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
128
     supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
129
     \f$-sup(u)\f$ demand.
130
     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}^+_0\f$
131
     solution of the following problem.
132 132

	
133
     \note A special case of this problem is when
134
     \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
135
     will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
136
     Thus a feasible solution for the
137
     \ref min_cost_flow "minimum cost flow" problem can be calculated
138
     in this way.
133
     \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
134
     \geq sup(u) \quad \forall u\in V, \f]
135
     \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
136
     
137
     The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
138
     zero or negative in order to have a feasible solution (since the sum
139
     of the expressions on the left-hand side of the inequalities is zero).
140
     It means that the total demand must be greater or equal to the total
141
     supply and all the supplies have to be carried out from the supply nodes,
142
     but there could be demands that are not satisfied.
143
     If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
144
     constraints have to be satisfied with equality, i.e. all demands
145
     have to be satisfied and all supplies have to be used.
146
     
147
     If you need the opposite inequalities in the supply/demand constraints
148
     (i.e. the total demand is less than the total supply and all the demands
149
     have to be satisfied while there could be supplies that are not used),
150
     then you could easily transform the problem to the above form by reversing
151
     the direction of the arcs and taking the negative of the supply values
152
     (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
153

	
154
     Note that this algorithm also provides a feasible solution for the
155
     \ref min_cost_flow "minimum cost flow problem".
139 156

	
140 157
     \tparam GR The type of the digraph the algorithm runs on.
141
     \tparam LM The type of the lower bound capacity map. The default
158
     \tparam LM The type of the lower bound map. The default
142 159
     map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
143
     \tparam UM The type of the upper bound capacity map. The default
144
     map type is \c LM.
145
     \tparam DM The type of the map that stores the lower bound
146
     for the supply of the nodes. The default map type is
160
     \tparam UM The type of the upper bound (capacity) map.
161
     The default map type is \c LM.
162
     \tparam SM The type of the supply map. The default map type is
147 163
     \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
148 164
  */
149 165
#ifdef DOXYGEN
150 166
template< typename GR,
151 167
          typename LM,
152 168
          typename UM,
153
          typename DM,
169
          typename SM,
154 170
          typename TR >
155 171
#else
156 172
template< typename GR,
157 173
          typename LM = typename GR::template ArcMap<int>,
158 174
          typename UM = LM,
159
          typename DM = typename GR::template NodeMap<typename UM::Value>,
160
          typename TR = CirculationDefaultTraits<GR, LM, UM, DM> >
175
          typename SM = typename GR::template NodeMap<typename UM::Value>,
176
          typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
161 177
#endif
162 178
  class Circulation {
163 179
  public:
164 180

	
165 181
    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
166 182
    typedef TR Traits;
167 183
    ///The type of the digraph the algorithm runs on.
168 184
    typedef typename Traits::Digraph Digraph;
169 185
    ///The type of the flow values.
170
    typedef typename Traits::Value Value;
186
    typedef typename Traits::Flow Flow;
171 187

	
172
    /// The type of the lower bound capacity map.
173
    typedef typename Traits::LCapMap LCapMap;
174
    /// The type of the upper bound capacity map.
175
    typedef typename Traits::UCapMap UCapMap;
176
    /// \brief The type of the map that stores the lower bound for
177
    /// the supply of the nodes.
178
    typedef typename Traits::DeltaMap DeltaMap;
188
    ///The type of the lower bound map.
189
    typedef typename Traits::LowerMap LowerMap;
190
    ///The type of the upper bound (capacity) map.
191
    typedef typename Traits::UpperMap UpperMap;
192
    ///The type of the supply map.
193
    typedef typename Traits::SupplyMap SupplyMap;
179 194
    ///The type of the flow map.
180 195
    typedef typename Traits::FlowMap FlowMap;
181 196

	
182 197
    ///The type of the elevator.
183 198
    typedef typename Traits::Elevator Elevator;
184 199
    ///The type of the tolerance.
185 200
    typedef typename Traits::Tolerance Tolerance;
186 201

	
187 202
  private:
188 203

	
189 204
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
190 205

	
191 206
    const Digraph &_g;
192 207
    int _node_num;
193 208

	
194
    const LCapMap *_lo;
195
    const UCapMap *_up;
196
    const DeltaMap *_delta;
209
    const LowerMap *_lo;
210
    const UpperMap *_up;
211
    const SupplyMap *_supply;
197 212

	
198 213
    FlowMap *_flow;
199 214
    bool _local_flow;
200 215

	
201 216
    Elevator* _level;
202 217
    bool _local_level;
203 218

	
204
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
219
    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
205 220
    ExcessMap* _excess;
206 221

	
207 222
    Tolerance _tol;
208 223
    int _el;
209 224

	
210 225
  public:
211 226

	
212 227
    typedef Circulation Create;
213 228

	
214 229
    ///\name Named Template Parameters
215 230

	
216 231
    ///@{
217 232

	
218 233
    template <typename T>
219 234
    struct SetFlowMapTraits : public Traits {
220 235
      typedef T FlowMap;
221 236
      static FlowMap *createFlowMap(const Digraph&) {
222 237
        LEMON_ASSERT(false, "FlowMap is not initialized");
223 238
        return 0; // ignore warnings
224 239
      }
225 240
    };
226 241

	
227 242
    /// \brief \ref named-templ-param "Named parameter" for setting
228 243
    /// FlowMap type
229 244
    ///
230 245
    /// \ref named-templ-param "Named parameter" for setting FlowMap
231 246
    /// type.
232 247
    template <typename T>
233 248
    struct SetFlowMap
234
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
249
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
235 250
                           SetFlowMapTraits<T> > {
236
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
251
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
237 252
                          SetFlowMapTraits<T> > Create;
238 253
    };
239 254

	
240 255
    template <typename T>
241 256
    struct SetElevatorTraits : public Traits {
242 257
      typedef T Elevator;
243 258
      static Elevator *createElevator(const Digraph&, int) {
244 259
        LEMON_ASSERT(false, "Elevator is not initialized");
245 260
        return 0; // ignore warnings
246 261
      }
247 262
    };
248 263

	
249 264
    /// \brief \ref named-templ-param "Named parameter" for setting
250 265
    /// Elevator type
251 266
    ///
252 267
    /// \ref named-templ-param "Named parameter" for setting Elevator
253 268
    /// type. If this named parameter is used, then an external
254 269
    /// elevator object must be passed to the algorithm using the
255 270
    /// \ref elevator(Elevator&) "elevator()" function before calling
256 271
    /// \ref run() or \ref init().
257 272
    /// \sa SetStandardElevator
258 273
    template <typename T>
259 274
    struct SetElevator
260
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
275
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
261 276
                           SetElevatorTraits<T> > {
262
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
277
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
263 278
                          SetElevatorTraits<T> > Create;
264 279
    };
265 280

	
266 281
    template <typename T>
267 282
    struct SetStandardElevatorTraits : public Traits {
268 283
      typedef T Elevator;
269 284
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
270 285
        return new Elevator(digraph, max_level);
271 286
      }
272 287
    };
273 288

	
274 289
    /// \brief \ref named-templ-param "Named parameter" for setting
275 290
    /// Elevator type with automatic allocation
276 291
    ///
277 292
    /// \ref named-templ-param "Named parameter" for setting Elevator
278 293
    /// type with automatic allocation.
279 294
    /// The Elevator should have standard constructor interface to be
280 295
    /// able to automatically created by the algorithm (i.e. the
281 296
    /// digraph and the maximum level should be passed to it).
282 297
    /// However an external elevator object could also be passed to the
283 298
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
284 299
    /// before calling \ref run() or \ref init().
285 300
    /// \sa SetElevator
286 301
    template <typename T>
287 302
    struct SetStandardElevator
288
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
303
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
289 304
                       SetStandardElevatorTraits<T> > {
290
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
305
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
291 306
                      SetStandardElevatorTraits<T> > Create;
292 307
    };
293 308

	
294 309
    /// @}
295 310

	
296 311
  protected:
297 312

	
298 313
    Circulation() {}
299 314

	
300 315
  public:
301 316

	
302
    /// The constructor of the class.
317
    /// Constructor.
303 318

	
304 319
    /// The constructor of the class.
305
    /// \param g The digraph the algorithm runs on.
306
    /// \param lo The lower bound capacity of the arcs.
307
    /// \param up The upper bound capacity of the arcs.
308
    /// \param delta The lower bound for the supply of the nodes.
309
    Circulation(const Digraph &g,const LCapMap &lo,
310
                const UCapMap &up,const DeltaMap &delta)
311
      : _g(g), _node_num(),
312
        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
313
        _level(0), _local_level(false), _excess(0), _el() {}
320
    ///
321
    /// \param graph The digraph the algorithm runs on.
322
    /// \param lower The lower bounds for the flow values on the arcs.
323
    /// \param upper The upper bounds (capacities) for the flow values 
324
    /// on the arcs.
325
    /// \param supply The signed supply values of the nodes.
326
    Circulation(const Digraph &graph, const LowerMap &lower,
327
                const UpperMap &upper, const SupplyMap &supply)
328
      : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
329
        _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
330
        _excess(NULL) {}
314 331

	
315 332
    /// Destructor.
316 333
    ~Circulation() {
317 334
      destroyStructures();
318 335
    }
319 336

	
320 337

	
321 338
  private:
322 339

	
323 340
    void createStructures() {
324 341
      _node_num = _el = countNodes(_g);
325 342

	
326 343
      if (!_flow) {
327 344
        _flow = Traits::createFlowMap(_g);
328 345
        _local_flow = true;
329 346
      }
330 347
      if (!_level) {
331 348
        _level = Traits::createElevator(_g, _node_num);
332 349
        _local_level = true;
333 350
      }
334 351
      if (!_excess) {
335 352
        _excess = new ExcessMap(_g);
336 353
      }
337 354
    }
338 355

	
339 356
    void destroyStructures() {
340 357
      if (_local_flow) {
341 358
        delete _flow;
342 359
      }
343 360
      if (_local_level) {
344 361
        delete _level;
345 362
      }
346 363
      if (_excess) {
347 364
        delete _excess;
348 365
      }
349 366
    }
350 367

	
351 368
  public:
352 369

	
353
    /// Sets the lower bound capacity map.
370
    /// Sets the lower bound map.
354 371

	
355
    /// Sets the lower bound capacity map.
372
    /// Sets the lower bound map.
356 373
    /// \return <tt>(*this)</tt>
357
    Circulation& lowerCapMap(const LCapMap& map) {
374
    Circulation& lowerMap(const LowerMap& map) {
358 375
      _lo = &map;
359 376
      return *this;
360 377
    }
361 378

	
362
    /// Sets the upper bound capacity map.
379
    /// Sets the upper bound (capacity) map.
363 380

	
364
    /// Sets the upper bound capacity map.
381
    /// Sets the upper bound (capacity) map.
365 382
    /// \return <tt>(*this)</tt>
366
    Circulation& upperCapMap(const LCapMap& map) {
383
    Circulation& upperMap(const LowerMap& map) {
367 384
      _up = &map;
368 385
      return *this;
369 386
    }
370 387

	
371
    /// Sets the lower bound map for the supply of the nodes.
388
    /// Sets the supply map.
372 389

	
373
    /// Sets the lower bound map for the supply of the nodes.
390
    /// Sets the supply map.
374 391
    /// \return <tt>(*this)</tt>
375
    Circulation& deltaMap(const DeltaMap& map) {
376
      _delta = &map;
392
    Circulation& supplyMap(const SupplyMap& map) {
393
      _supply = &map;
377 394
      return *this;
378 395
    }
379 396

	
380 397
    /// \brief Sets the flow map.
381 398
    ///
382 399
    /// Sets the flow map.
383 400
    /// If you don't use this function before calling \ref run() or
384 401
    /// \ref init(), an instance will be allocated automatically.
385 402
    /// The destructor deallocates this automatically allocated map,
386 403
    /// of course.
387 404
    /// \return <tt>(*this)</tt>
388 405
    Circulation& flowMap(FlowMap& map) {
389 406
      if (_local_flow) {
390 407
        delete _flow;
391 408
        _local_flow = false;
392 409
      }
393 410
      _flow = &map;
394 411
      return *this;
395 412
    }
396 413

	
397 414
    /// \brief Sets the elevator used by algorithm.
398 415
    ///
399 416
    /// Sets the elevator used by algorithm.
400 417
    /// If you don't use this function before calling \ref run() or
401 418
    /// \ref init(), an instance will be allocated automatically.
402 419
    /// The destructor deallocates this automatically allocated elevator,
403 420
    /// of course.
404 421
    /// \return <tt>(*this)</tt>
405 422
    Circulation& elevator(Elevator& elevator) {
406 423
      if (_local_level) {
407 424
        delete _level;
408 425
        _local_level = false;
409 426
      }
410 427
      _level = &elevator;
411 428
      return *this;
412 429
    }
413 430

	
414 431
    /// \brief Returns a const reference to the elevator.
415 432
    ///
416 433
    /// Returns a const reference to the elevator.
417 434
    ///
418 435
    /// \pre Either \ref run() or \ref init() must be called before
419 436
    /// using this function.
420 437
    const Elevator& elevator() const {
421 438
      return *_level;
422 439
    }
423 440

	
424 441
    /// \brief Sets the tolerance used by algorithm.
425 442
    ///
426 443
    /// Sets the tolerance used by algorithm.
427 444
    Circulation& tolerance(const Tolerance& tolerance) const {
428 445
      _tol = tolerance;
429 446
      return *this;
430 447
    }
431 448

	
432 449
    /// \brief Returns a const reference to the tolerance.
433 450
    ///
434 451
    /// Returns a const reference to the tolerance.
435 452
    const Tolerance& tolerance() const {
436 453
      return tolerance;
437 454
    }
438 455

	
439 456
    /// \name Execution Control
440 457
    /// The simplest way to execute the algorithm is to call \ref run().\n
441 458
    /// If you need more control on the initial solution or the execution,
442 459
    /// first you have to call one of the \ref init() functions, then
443 460
    /// the \ref start() function.
444 461

	
445 462
    ///@{
446 463

	
447 464
    /// Initializes the internal data structures.
448 465

	
449 466
    /// Initializes the internal data structures and sets all flow values
450 467
    /// to the lower bound.
451 468
    void init()
452 469
    {
453 470
      createStructures();
454 471

	
455 472
      for(NodeIt n(_g);n!=INVALID;++n) {
456
        (*_excess)[n] = (*_delta)[n];
473
        (*_excess)[n] = (*_supply)[n];
457 474
      }
458 475

	
459 476
      for (ArcIt e(_g);e!=INVALID;++e) {
460 477
        _flow->set(e, (*_lo)[e]);
461 478
        (*_excess)[_g.target(e)] += (*_flow)[e];
462 479
        (*_excess)[_g.source(e)] -= (*_flow)[e];
463 480
      }
464 481

	
465 482
      // global relabeling tested, but in general case it provides
466 483
      // worse performance for random digraphs
467 484
      _level->initStart();
468 485
      for(NodeIt n(_g);n!=INVALID;++n)
469 486
        _level->initAddItem(n);
470 487
      _level->initFinish();
471 488
      for(NodeIt n(_g);n!=INVALID;++n)
472 489
        if(_tol.positive((*_excess)[n]))
473 490
          _level->activate(n);
474 491
    }
475 492

	
476 493
    /// Initializes the internal data structures using a greedy approach.
477 494

	
478 495
    /// Initializes the internal data structures using a greedy approach
479 496
    /// to construct the initial solution.
480 497
    void greedyInit()
481 498
    {
482 499
      createStructures();
483 500

	
484 501
      for(NodeIt n(_g);n!=INVALID;++n) {
485
        (*_excess)[n] = (*_delta)[n];
502
        (*_excess)[n] = (*_supply)[n];
486 503
      }
487 504

	
488 505
      for (ArcIt e(_g);e!=INVALID;++e) {
489 506
        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
490 507
          _flow->set(e, (*_up)[e]);
491 508
          (*_excess)[_g.target(e)] += (*_up)[e];
492 509
          (*_excess)[_g.source(e)] -= (*_up)[e];
493 510
        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
494 511
          _flow->set(e, (*_lo)[e]);
495 512
          (*_excess)[_g.target(e)] += (*_lo)[e];
496 513
          (*_excess)[_g.source(e)] -= (*_lo)[e];
497 514
        } else {
498
          Value fc = -(*_excess)[_g.target(e)];
515
          Flow fc = -(*_excess)[_g.target(e)];
499 516
          _flow->set(e, fc);
500 517
          (*_excess)[_g.target(e)] = 0;
501 518
          (*_excess)[_g.source(e)] -= fc;
502 519
        }
503 520
      }
504 521

	
505 522
      _level->initStart();
506 523
      for(NodeIt n(_g);n!=INVALID;++n)
507 524
        _level->initAddItem(n);
508 525
      _level->initFinish();
509 526
      for(NodeIt n(_g);n!=INVALID;++n)
510 527
        if(_tol.positive((*_excess)[n]))
511 528
          _level->activate(n);
512 529
    }
513 530

	
514 531
    ///Executes the algorithm
515 532

	
516 533
    ///This function executes the algorithm.
517 534
    ///
518 535
    ///\return \c true if a feasible circulation is found.
519 536
    ///
520 537
    ///\sa barrier()
521 538
    ///\sa barrierMap()
522 539
    bool start()
523 540
    {
524 541

	
525 542
      Node act;
526 543
      Node bact=INVALID;
527 544
      Node last_activated=INVALID;
528 545
      while((act=_level->highestActive())!=INVALID) {
529 546
        int actlevel=(*_level)[act];
530 547
        int mlevel=_node_num;
531
        Value exc=(*_excess)[act];
548
        Flow exc=(*_excess)[act];
532 549

	
533 550
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
534 551
          Node v = _g.target(e);
535
          Value fc=(*_up)[e]-(*_flow)[e];
552
          Flow fc=(*_up)[e]-(*_flow)[e];
536 553
          if(!_tol.positive(fc)) continue;
537 554
          if((*_level)[v]<actlevel) {
538 555
            if(!_tol.less(fc, exc)) {
539 556
              _flow->set(e, (*_flow)[e] + exc);
540 557
              (*_excess)[v] += exc;
541 558
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
542 559
                _level->activate(v);
543 560
              (*_excess)[act] = 0;
544 561
              _level->deactivate(act);
545 562
              goto next_l;
546 563
            }
547 564
            else {
548 565
              _flow->set(e, (*_up)[e]);
549 566
              (*_excess)[v] += fc;
550 567
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
551 568
                _level->activate(v);
552 569
              exc-=fc;
553 570
            }
554 571
          }
555 572
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
556 573
        }
557 574
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
558 575
          Node v = _g.source(e);
559
          Value fc=(*_flow)[e]-(*_lo)[e];
576
          Flow fc=(*_flow)[e]-(*_lo)[e];
560 577
          if(!_tol.positive(fc)) continue;
561 578
          if((*_level)[v]<actlevel) {
562 579
            if(!_tol.less(fc, exc)) {
563 580
              _flow->set(e, (*_flow)[e] - exc);
564 581
              (*_excess)[v] += exc;
565 582
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
566 583
                _level->activate(v);
567 584
              (*_excess)[act] = 0;
568 585
              _level->deactivate(act);
569 586
              goto next_l;
570 587
            }
571 588
            else {
572 589
              _flow->set(e, (*_lo)[e]);
573 590
              (*_excess)[v] += fc;
574 591
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
575 592
                _level->activate(v);
576 593
              exc-=fc;
577 594
            }
578 595
          }
579 596
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
580 597
        }
581 598

	
582 599
        (*_excess)[act] = exc;
583 600
        if(!_tol.positive(exc)) _level->deactivate(act);
584 601
        else if(mlevel==_node_num) {
585 602
          _level->liftHighestActiveToTop();
586 603
          _el = _node_num;
587 604
          return false;
588 605
        }
589 606
        else {
590 607
          _level->liftHighestActive(mlevel+1);
591 608
          if(_level->onLevel(actlevel)==0) {
592 609
            _el = actlevel;
593 610
            return false;
594 611
          }
595 612
        }
596 613
      next_l:
597 614
        ;
598 615
      }
599 616
      return true;
600 617
    }
601 618

	
602 619
    /// Runs the algorithm.
603 620

	
604 621
    /// This function runs the algorithm.
605 622
    ///
606 623
    /// \return \c true if a feasible circulation is found.
607 624
    ///
608 625
    /// \note Apart from the return value, c.run() is just a shortcut of
609 626
    /// the following code.
610 627
    /// \code
611 628
    ///   c.greedyInit();
612 629
    ///   c.start();
613 630
    /// \endcode
614 631
    bool run() {
615 632
      greedyInit();
616 633
      return start();
617 634
    }
618 635

	
619 636
    /// @}
620 637

	
621 638
    /// \name Query Functions
622 639
    /// The results of the circulation algorithm can be obtained using
623 640
    /// these functions.\n
624 641
    /// Either \ref run() or \ref start() should be called before
625 642
    /// using them.
626 643

	
627 644
    ///@{
628 645

	
629 646
    /// \brief Returns the flow on the given arc.
630 647
    ///
631 648
    /// Returns the flow on the given arc.
632 649
    ///
633 650
    /// \pre Either \ref run() or \ref init() must be called before
634 651
    /// using this function.
635
    Value flow(const Arc& arc) const {
652
    Flow flow(const Arc& arc) const {
636 653
      return (*_flow)[arc];
637 654
    }
638 655

	
639 656
    /// \brief Returns a const reference to the flow map.
640 657
    ///
641 658
    /// Returns a const reference to the arc map storing the found flow.
642 659
    ///
643 660
    /// \pre Either \ref run() or \ref init() must be called before
644 661
    /// using this function.
645 662
    const FlowMap& flowMap() const {
646 663
      return *_flow;
647 664
    }
648 665

	
649 666
    /**
650 667
       \brief Returns \c true if the given node is in a barrier.
651 668

	
652 669
       Barrier is a set \e B of nodes for which
653 670

	
654
       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
655
           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
671
       \f[ \sum_{uv\in A: u\in B} upper(uv) -
672
           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
656 673

	
657 674
       holds. The existence of a set with this property prooves that a
658 675
       feasible circualtion cannot exist.
659 676

	
660 677
       This function returns \c true if the given node is in the found
661 678
       barrier. If a feasible circulation is found, the function
662 679
       gives back \c false for every node.
663 680

	
664 681
       \pre Either \ref run() or \ref init() must be called before
665 682
       using this function.
666 683

	
667 684
       \sa barrierMap()
668 685
       \sa checkBarrier()
669 686
    */
670 687
    bool barrier(const Node& node) const
671 688
    {
672 689
      return (*_level)[node] >= _el;
673 690
    }
674 691

	
675 692
    /// \brief Gives back a barrier.
676 693
    ///
677 694
    /// This function sets \c bar to the characteristic vector of the
678 695
    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
679 696
    /// node map with \c bool (or convertible) value type.
680 697
    ///
681 698
    /// If a feasible circulation is found, the function gives back an
682 699
    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
683 700
    ///
684 701
    /// \note This function calls \ref barrier() for each node,
685 702
    /// so it runs in O(n) time.
686 703
    ///
687 704
    /// \pre Either \ref run() or \ref init() must be called before
688 705
    /// using this function.
689 706
    ///
690 707
    /// \sa barrier()
691 708
    /// \sa checkBarrier()
692 709
    template<class BarrierMap>
693 710
    void barrierMap(BarrierMap &bar) const
694 711
    {
695 712
      for(NodeIt n(_g);n!=INVALID;++n)
696 713
        bar.set(n, (*_level)[n] >= _el);
697 714
    }
698 715

	
699 716
    /// @}
700 717

	
701 718
    /// \name Checker Functions
702 719
    /// The feasibility of the results can be checked using
703 720
    /// these functions.\n
704 721
    /// Either \ref run() or \ref start() should be called before
705 722
    /// using them.
706 723

	
707 724
    ///@{
708 725

	
709 726
    ///Check if the found flow is a feasible circulation
710 727

	
711 728
    ///Check if the found flow is a feasible circulation,
712 729
    ///
713 730
    bool checkFlow() const {
714 731
      for(ArcIt e(_g);e!=INVALID;++e)
715 732
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
716 733
      for(NodeIt n(_g);n!=INVALID;++n)
717 734
        {
718
          Value dif=-(*_delta)[n];
735
          Flow dif=-(*_supply)[n];
719 736
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
720 737
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
721 738
          if(_tol.negative(dif)) return false;
722 739
        }
723 740
      return true;
724 741
    }
725 742

	
726 743
    ///Check whether or not the last execution provides a barrier
727 744

	
728 745
    ///Check whether or not the last execution provides a barrier.
729 746
    ///\sa barrier()
730 747
    ///\sa barrierMap()
731 748
    bool checkBarrier() const
732 749
    {
733
      Value delta=0;
750
      Flow delta=0;
734 751
      for(NodeIt n(_g);n!=INVALID;++n)
735 752
        if(barrier(n))
736
          delta-=(*_delta)[n];
753
          delta-=(*_supply)[n];
737 754
      for(ArcIt e(_g);e!=INVALID;++e)
738 755
        {
739 756
          Node s=_g.source(e);
740 757
          Node t=_g.target(e);
741 758
          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
742 759
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
743 760
        }
744 761
      return _tol.negative(delta);
745 762
    }
746 763

	
747 764
    /// @}
748 765

	
749 766
  };
750 767

	
751 768
}
752 769

	
753 770
#endif
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_PREFLOW_H
20 20
#define LEMON_PREFLOW_H
21 21

	
22 22
#include <lemon/tolerance.h>
23 23
#include <lemon/elevator.h>
24 24

	
25 25
/// \file
26 26
/// \ingroup max_flow
27 27
/// \brief Implementation of the preflow algorithm.
28 28

	
29 29
namespace lemon {
30 30

	
31 31
  /// \brief Default traits class of Preflow class.
32 32
  ///
33 33
  /// Default traits class of Preflow class.
34 34
  /// \tparam GR Digraph type.
35 35
  /// \tparam CAP Capacity map type.
36 36
  template <typename GR, typename CAP>
37 37
  struct PreflowDefaultTraits {
38 38

	
39 39
    /// \brief The type of the digraph the algorithm runs on.
40 40
    typedef GR Digraph;
41 41

	
42 42
    /// \brief The type of the map that stores the arc capacities.
43 43
    ///
44 44
    /// The type of the map that stores the arc capacities.
45 45
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
46 46
    typedef CAP CapacityMap;
47 47

	
48 48
    /// \brief The type of the flow values.
49
    typedef typename CapacityMap::Value Value;
49
    typedef typename CapacityMap::Value Flow;
50 50

	
51 51
    /// \brief The type of the map that stores the flow values.
52 52
    ///
53 53
    /// The type of the map that stores the flow values.
54 54
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
55
    typedef typename Digraph::template ArcMap<Value> FlowMap;
55
    typedef typename Digraph::template ArcMap<Flow> FlowMap;
56 56

	
57 57
    /// \brief Instantiates a FlowMap.
58 58
    ///
59 59
    /// This function instantiates a \ref FlowMap.
60
    /// \param digraph The digraph, to which we would like to define
60
    /// \param digraph The digraph for which we would like to define
61 61
    /// the flow map.
62 62
    static FlowMap* createFlowMap(const Digraph& digraph) {
63 63
      return new FlowMap(digraph);
64 64
    }
65 65

	
66 66
    /// \brief The elevator type used by Preflow algorithm.
67 67
    ///
68 68
    /// The elevator type used by Preflow algorithm.
69 69
    ///
70 70
    /// \sa Elevator
71 71
    /// \sa LinkedElevator
72 72
    typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
73 73

	
74 74
    /// \brief Instantiates an Elevator.
75 75
    ///
76 76
    /// This function instantiates an \ref Elevator.
77
    /// \param digraph The digraph, to which we would like to define
77
    /// \param digraph The digraph for which we would like to define
78 78
    /// the elevator.
79 79
    /// \param max_level The maximum level of the elevator.
80 80
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
81 81
      return new Elevator(digraph, max_level);
82 82
    }
83 83

	
84 84
    /// \brief The tolerance used by the algorithm
85 85
    ///
86 86
    /// The tolerance used by the algorithm to handle inexact computation.
87
    typedef lemon::Tolerance<Value> Tolerance;
87
    typedef lemon::Tolerance<Flow> Tolerance;
88 88

	
89 89
  };
90 90

	
91 91

	
92 92
  /// \ingroup max_flow
93 93
  ///
94 94
  /// \brief %Preflow algorithm class.
95 95
  ///
96 96
  /// This class provides an implementation of Goldberg-Tarjan's \e preflow
97 97
  /// \e push-relabel algorithm producing a \ref max_flow
98 98
  /// "flow of maximum value" in a digraph.
99 99
  /// The preflow algorithms are the fastest known maximum
100 100
  /// flow algorithms. The current implementation use a mixture of the
101 101
  /// \e "highest label" and the \e "bound decrease" heuristics.
102 102
  /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
103 103
  ///
104 104
  /// The algorithm consists of two phases. After the first phase
105 105
  /// the maximum flow value and the minimum cut is obtained. The
106 106
  /// second phase constructs a feasible maximum flow on each arc.
107 107
  ///
108 108
  /// \tparam GR The type of the digraph the algorithm runs on.
109 109
  /// \tparam CAP The type of the capacity map. The default map
110 110
  /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
111 111
#ifdef DOXYGEN
112 112
  template <typename GR, typename CAP, typename TR>
113 113
#else
114 114
  template <typename GR,
115 115
            typename CAP = typename GR::template ArcMap<int>,
116 116
            typename TR = PreflowDefaultTraits<GR, CAP> >
117 117
#endif
118 118
  class Preflow {
119 119
  public:
120 120

	
121 121
    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
122 122
    typedef TR Traits;
123 123
    ///The type of the digraph the algorithm runs on.
124 124
    typedef typename Traits::Digraph Digraph;
125 125
    ///The type of the capacity map.
126 126
    typedef typename Traits::CapacityMap CapacityMap;
127 127
    ///The type of the flow values.
128
    typedef typename Traits::Value Value;
128
    typedef typename Traits::Flow Flow;
129 129

	
130 130
    ///The type of the flow map.
131 131
    typedef typename Traits::FlowMap FlowMap;
132 132
    ///The type of the elevator.
133 133
    typedef typename Traits::Elevator Elevator;
134 134
    ///The type of the tolerance.
135 135
    typedef typename Traits::Tolerance Tolerance;
136 136

	
137 137
  private:
138 138

	
139 139
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
140 140

	
141 141
    const Digraph& _graph;
142 142
    const CapacityMap* _capacity;
143 143

	
144 144
    int _node_num;
145 145

	
146 146
    Node _source, _target;
147 147

	
148 148
    FlowMap* _flow;
149 149
    bool _local_flow;
150 150

	
151 151
    Elevator* _level;
152 152
    bool _local_level;
153 153

	
154
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
154
    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
155 155
    ExcessMap* _excess;
156 156

	
157 157
    Tolerance _tolerance;
158 158

	
159 159
    bool _phase;
160 160

	
161 161

	
162 162
    void createStructures() {
163 163
      _node_num = countNodes(_graph);
164 164

	
165 165
      if (!_flow) {
166 166
        _flow = Traits::createFlowMap(_graph);
167 167
        _local_flow = true;
168 168
      }
169 169
      if (!_level) {
170 170
        _level = Traits::createElevator(_graph, _node_num);
171 171
        _local_level = true;
172 172
      }
173 173
      if (!_excess) {
174 174
        _excess = new ExcessMap(_graph);
175 175
      }
176 176
    }
177 177

	
178 178
    void destroyStructures() {
179 179
      if (_local_flow) {
180 180
        delete _flow;
181 181
      }
182 182
      if (_local_level) {
183 183
        delete _level;
184 184
      }
185 185
      if (_excess) {
186 186
        delete _excess;
187 187
      }
188 188
    }
189 189

	
190 190
  public:
191 191

	
192 192
    typedef Preflow Create;
193 193

	
194 194
    ///\name Named Template Parameters
195 195

	
196 196
    ///@{
197 197

	
198 198
    template <typename T>
199 199
    struct SetFlowMapTraits : public Traits {
200 200
      typedef T FlowMap;
201 201
      static FlowMap *createFlowMap(const Digraph&) {
202 202
        LEMON_ASSERT(false, "FlowMap is not initialized");
203 203
        return 0; // ignore warnings
204 204
      }
205 205
    };
206 206

	
207 207
    /// \brief \ref named-templ-param "Named parameter" for setting
208 208
    /// FlowMap type
209 209
    ///
210 210
    /// \ref named-templ-param "Named parameter" for setting FlowMap
211 211
    /// type.
212 212
    template <typename T>
213 213
    struct SetFlowMap
214 214
      : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > {
215 215
      typedef Preflow<Digraph, CapacityMap,
216 216
                      SetFlowMapTraits<T> > Create;
217 217
    };
218 218

	
219 219
    template <typename T>
220 220
    struct SetElevatorTraits : public Traits {
221 221
      typedef T Elevator;
222 222
      static Elevator *createElevator(const Digraph&, int) {
223 223
        LEMON_ASSERT(false, "Elevator is not initialized");
224 224
        return 0; // ignore warnings
225 225
      }
226 226
    };
227 227

	
228 228
    /// \brief \ref named-templ-param "Named parameter" for setting
229 229
    /// Elevator type
230 230
    ///
231 231
    /// \ref named-templ-param "Named parameter" for setting Elevator
232 232
    /// type. If this named parameter is used, then an external
233 233
    /// elevator object must be passed to the algorithm using the
234 234
    /// \ref elevator(Elevator&) "elevator()" function before calling
235 235
    /// \ref run() or \ref init().
236 236
    /// \sa SetStandardElevator
237 237
    template <typename T>
238 238
    struct SetElevator
239 239
      : public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > {
240 240
      typedef Preflow<Digraph, CapacityMap,
241 241
                      SetElevatorTraits<T> > Create;
242 242
    };
243 243

	
244 244
    template <typename T>
245 245
    struct SetStandardElevatorTraits : public Traits {
246 246
      typedef T Elevator;
247 247
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
248 248
        return new Elevator(digraph, max_level);
249 249
      }
250 250
    };
251 251

	
252 252
    /// \brief \ref named-templ-param "Named parameter" for setting
253 253
    /// Elevator type with automatic allocation
254 254
    ///
255 255
    /// \ref named-templ-param "Named parameter" for setting Elevator
256 256
    /// type with automatic allocation.
257 257
    /// The Elevator should have standard constructor interface to be
258 258
    /// able to automatically created by the algorithm (i.e. the
259 259
    /// digraph and the maximum level should be passed to it).
260 260
    /// However an external elevator object could also be passed to the
261 261
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
262 262
    /// before calling \ref run() or \ref init().
263 263
    /// \sa SetElevator
264 264
    template <typename T>
265 265
    struct SetStandardElevator
266 266
      : public Preflow<Digraph, CapacityMap,
267 267
                       SetStandardElevatorTraits<T> > {
268 268
      typedef Preflow<Digraph, CapacityMap,
269 269
                      SetStandardElevatorTraits<T> > Create;
270 270
    };
271 271

	
272 272
    /// @}
273 273

	
274 274
  protected:
275 275

	
276 276
    Preflow() {}
277 277

	
278 278
  public:
279 279

	
280 280

	
281 281
    /// \brief The constructor of the class.
282 282
    ///
283 283
    /// The constructor of the class.
284 284
    /// \param digraph The digraph the algorithm runs on.
285 285
    /// \param capacity The capacity of the arcs.
286 286
    /// \param source The source node.
287 287
    /// \param target The target node.
288 288
    Preflow(const Digraph& digraph, const CapacityMap& capacity,
289 289
            Node source, Node target)
290 290
      : _graph(digraph), _capacity(&capacity),
291 291
        _node_num(0), _source(source), _target(target),
292 292
        _flow(0), _local_flow(false),
293 293
        _level(0), _local_level(false),
294 294
        _excess(0), _tolerance(), _phase() {}
295 295

	
296 296
    /// \brief Destructor.
297 297
    ///
298 298
    /// Destructor.
299 299
    ~Preflow() {
300 300
      destroyStructures();
301 301
    }
302 302

	
303 303
    /// \brief Sets the capacity map.
304 304
    ///
305 305
    /// Sets the capacity map.
306 306
    /// \return <tt>(*this)</tt>
307 307
    Preflow& capacityMap(const CapacityMap& map) {
308 308
      _capacity = &map;
309 309
      return *this;
310 310
    }
311 311

	
312 312
    /// \brief Sets the flow map.
313 313
    ///
314 314
    /// Sets the flow map.
315 315
    /// If you don't use this function before calling \ref run() or
316 316
    /// \ref init(), an instance will be allocated automatically.
317 317
    /// The destructor deallocates this automatically allocated map,
318 318
    /// of course.
319 319
    /// \return <tt>(*this)</tt>
320 320
    Preflow& flowMap(FlowMap& map) {
321 321
      if (_local_flow) {
322 322
        delete _flow;
323 323
        _local_flow = false;
324 324
      }
325 325
      _flow = &map;
326 326
      return *this;
327 327
    }
328 328

	
329 329
    /// \brief Sets the source node.
330 330
    ///
331 331
    /// Sets the source node.
332 332
    /// \return <tt>(*this)</tt>
333 333
    Preflow& source(const Node& node) {
334 334
      _source = node;
335 335
      return *this;
336 336
    }
337 337

	
338 338
    /// \brief Sets the target node.
339 339
    ///
340 340
    /// Sets the target node.
341 341
    /// \return <tt>(*this)</tt>
342 342
    Preflow& target(const Node& node) {
343 343
      _target = node;
344 344
      return *this;
345 345
    }
346 346

	
347 347
    /// \brief Sets the elevator used by algorithm.
348 348
    ///
349 349
    /// Sets the elevator used by algorithm.
350 350
    /// If you don't use this function before calling \ref run() or
351 351
    /// \ref init(), an instance will be allocated automatically.
352 352
    /// The destructor deallocates this automatically allocated elevator,
353 353
    /// of course.
354 354
    /// \return <tt>(*this)</tt>
355 355
    Preflow& elevator(Elevator& elevator) {
356 356
      if (_local_level) {
357 357
        delete _level;
358 358
        _local_level = false;
359 359
      }
360 360
      _level = &elevator;
361 361
      return *this;
362 362
    }
363 363

	
364 364
    /// \brief Returns a const reference to the elevator.
365 365
    ///
366 366
    /// Returns a const reference to the elevator.
367 367
    ///
368 368
    /// \pre Either \ref run() or \ref init() must be called before
369 369
    /// using this function.
370 370
    const Elevator& elevator() const {
371 371
      return *_level;
372 372
    }
373 373

	
374 374
    /// \brief Sets the tolerance used by algorithm.
375 375
    ///
376 376
    /// Sets the tolerance used by algorithm.
377 377
    Preflow& tolerance(const Tolerance& tolerance) const {
378 378
      _tolerance = tolerance;
379 379
      return *this;
380 380
    }
381 381

	
382 382
    /// \brief Returns a const reference to the tolerance.
383 383
    ///
384 384
    /// Returns a const reference to the tolerance.
385 385
    const Tolerance& tolerance() const {
386 386
      return tolerance;
387 387
    }
388 388

	
389 389
    /// \name Execution Control
390 390
    /// The simplest way to execute the preflow algorithm is to use
391 391
    /// \ref run() or \ref runMinCut().\n
392 392
    /// If you need more control on the initial solution or the execution,
393 393
    /// first you have to call one of the \ref init() functions, then
394 394
    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
395 395

	
396 396
    ///@{
397 397

	
398 398
    /// \brief Initializes the internal data structures.
399 399
    ///
400 400
    /// Initializes the internal data structures and sets the initial
401 401
    /// flow to zero on each arc.
402 402
    void init() {
403 403
      createStructures();
404 404

	
405 405
      _phase = true;
406 406
      for (NodeIt n(_graph); n != INVALID; ++n) {
407 407
        (*_excess)[n] = 0;
408 408
      }
409 409

	
410 410
      for (ArcIt e(_graph); e != INVALID; ++e) {
411 411
        _flow->set(e, 0);
412 412
      }
413 413

	
414 414
      typename Digraph::template NodeMap<bool> reached(_graph, false);
415 415

	
416 416
      _level->initStart();
417 417
      _level->initAddItem(_target);
418 418

	
419 419
      std::vector<Node> queue;
420 420
      reached[_source] = true;
421 421

	
422 422
      queue.push_back(_target);
423 423
      reached[_target] = true;
424 424
      while (!queue.empty()) {
425 425
        _level->initNewLevel();
426 426
        std::vector<Node> nqueue;
427 427
        for (int i = 0; i < int(queue.size()); ++i) {
428 428
          Node n = queue[i];
429 429
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
430 430
            Node u = _graph.source(e);
431 431
            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
432 432
              reached[u] = true;
433 433
              _level->initAddItem(u);
434 434
              nqueue.push_back(u);
435 435
            }
436 436
          }
437 437
        }
438 438
        queue.swap(nqueue);
439 439
      }
440 440
      _level->initFinish();
441 441

	
442 442
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
443 443
        if (_tolerance.positive((*_capacity)[e])) {
444 444
          Node u = _graph.target(e);
445 445
          if ((*_level)[u] == _level->maxLevel()) continue;
446 446
          _flow->set(e, (*_capacity)[e]);
447 447
          (*_excess)[u] += (*_capacity)[e];
448 448
          if (u != _target && !_level->active(u)) {
449 449
            _level->activate(u);
450 450
          }
451 451
        }
452 452
      }
453 453
    }
454 454

	
455 455
    /// \brief Initializes the internal data structures using the
456 456
    /// given flow map.
457 457
    ///
458 458
    /// Initializes the internal data structures and sets the initial
459 459
    /// flow to the given \c flowMap. The \c flowMap should contain a
460 460
    /// flow or at least a preflow, i.e. at each node excluding the
461 461
    /// source node the incoming flow should greater or equal to the
462 462
    /// outgoing flow.
463 463
    /// \return \c false if the given \c flowMap is not a preflow.
464 464
    template <typename FlowMap>
465 465
    bool init(const FlowMap& flowMap) {
466 466
      createStructures();
467 467

	
468 468
      for (ArcIt e(_graph); e != INVALID; ++e) {
469 469
        _flow->set(e, flowMap[e]);
470 470
      }
471 471

	
472 472
      for (NodeIt n(_graph); n != INVALID; ++n) {
473
        Value excess = 0;
473
        Flow excess = 0;
474 474
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
475 475
          excess += (*_flow)[e];
476 476
        }
477 477
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
478 478
          excess -= (*_flow)[e];
479 479
        }
480 480
        if (excess < 0 && n != _source) return false;
481 481
        (*_excess)[n] = excess;
482 482
      }
483 483

	
484 484
      typename Digraph::template NodeMap<bool> reached(_graph, false);
485 485

	
486 486
      _level->initStart();
487 487
      _level->initAddItem(_target);
488 488

	
489 489
      std::vector<Node> queue;
490 490
      reached[_source] = true;
491 491

	
492 492
      queue.push_back(_target);
493 493
      reached[_target] = true;
494 494
      while (!queue.empty()) {
495 495
        _level->initNewLevel();
496 496
        std::vector<Node> nqueue;
497 497
        for (int i = 0; i < int(queue.size()); ++i) {
498 498
          Node n = queue[i];
499 499
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
500 500
            Node u = _graph.source(e);
501 501
            if (!reached[u] &&
502 502
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
503 503
              reached[u] = true;
504 504
              _level->initAddItem(u);
505 505
              nqueue.push_back(u);
506 506
            }
507 507
          }
508 508
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
509 509
            Node v = _graph.target(e);
510 510
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
511 511
              reached[v] = true;
512 512
              _level->initAddItem(v);
513 513
              nqueue.push_back(v);
514 514
            }
515 515
          }
516 516
        }
517 517
        queue.swap(nqueue);
518 518
      }
519 519
      _level->initFinish();
520 520

	
521 521
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
522
        Value rem = (*_capacity)[e] - (*_flow)[e];
522
        Flow rem = (*_capacity)[e] - (*_flow)[e];
523 523
        if (_tolerance.positive(rem)) {
524 524
          Node u = _graph.target(e);
525 525
          if ((*_level)[u] == _level->maxLevel()) continue;
526 526
          _flow->set(e, (*_capacity)[e]);
527 527
          (*_excess)[u] += rem;
528 528
          if (u != _target && !_level->active(u)) {
529 529
            _level->activate(u);
530 530
          }
531 531
        }
532 532
      }
533 533
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
534
        Value rem = (*_flow)[e];
534
        Flow rem = (*_flow)[e];
535 535
        if (_tolerance.positive(rem)) {
536 536
          Node v = _graph.source(e);
537 537
          if ((*_level)[v] == _level->maxLevel()) continue;
538 538
          _flow->set(e, 0);
539 539
          (*_excess)[v] += rem;
540 540
          if (v != _target && !_level->active(v)) {
541 541
            _level->activate(v);
542 542
          }
543 543
        }
544 544
      }
545 545
      return true;
546 546
    }
547 547

	
548 548
    /// \brief Starts the first phase of the preflow algorithm.
549 549
    ///
550 550
    /// The preflow algorithm consists of two phases, this method runs
551 551
    /// the first phase. After the first phase the maximum flow value
552 552
    /// and a minimum value cut can already be computed, although a
553 553
    /// maximum flow is not yet obtained. So after calling this method
554 554
    /// \ref flowValue() returns the value of a maximum flow and \ref
555 555
    /// minCut() returns a minimum cut.
556 556
    /// \pre One of the \ref init() functions must be called before
557 557
    /// using this function.
558 558
    void startFirstPhase() {
559 559
      _phase = true;
560 560

	
561 561
      Node n = _level->highestActive();
562 562
      int level = _level->highestActiveLevel();
563 563
      while (n != INVALID) {
564 564
        int num = _node_num;
565 565

	
566 566
        while (num > 0 && n != INVALID) {
567
          Value excess = (*_excess)[n];
567
          Flow excess = (*_excess)[n];
568 568
          int new_level = _level->maxLevel();
569 569

	
570 570
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
571
            Value rem = (*_capacity)[e] - (*_flow)[e];
571
            Flow rem = (*_capacity)[e] - (*_flow)[e];
572 572
            if (!_tolerance.positive(rem)) continue;
573 573
            Node v = _graph.target(e);
574 574
            if ((*_level)[v] < level) {
575 575
              if (!_level->active(v) && v != _target) {
576 576
                _level->activate(v);
577 577
              }
578 578
              if (!_tolerance.less(rem, excess)) {
579 579
                _flow->set(e, (*_flow)[e] + excess);
580 580
                (*_excess)[v] += excess;
581 581
                excess = 0;
582 582
                goto no_more_push_1;
583 583
              } else {
584 584
                excess -= rem;
585 585
                (*_excess)[v] += rem;
586 586
                _flow->set(e, (*_capacity)[e]);
587 587
              }
588 588
            } else if (new_level > (*_level)[v]) {
589 589
              new_level = (*_level)[v];
590 590
            }
591 591
          }
592 592

	
593 593
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
594
            Value rem = (*_flow)[e];
594
            Flow rem = (*_flow)[e];
595 595
            if (!_tolerance.positive(rem)) continue;
596 596
            Node v = _graph.source(e);
597 597
            if ((*_level)[v] < level) {
598 598
              if (!_level->active(v) && v != _target) {
599 599
                _level->activate(v);
600 600
              }
601 601
              if (!_tolerance.less(rem, excess)) {
602 602
                _flow->set(e, (*_flow)[e] - excess);
603 603
                (*_excess)[v] += excess;
604 604
                excess = 0;
605 605
                goto no_more_push_1;
606 606
              } else {
607 607
                excess -= rem;
608 608
                (*_excess)[v] += rem;
609 609
                _flow->set(e, 0);
610 610
              }
611 611
            } else if (new_level > (*_level)[v]) {
612 612
              new_level = (*_level)[v];
613 613
            }
614 614
          }
615 615

	
616 616
        no_more_push_1:
617 617

	
618 618
          (*_excess)[n] = excess;
619 619

	
620 620
          if (excess != 0) {
621 621
            if (new_level + 1 < _level->maxLevel()) {
622 622
              _level->liftHighestActive(new_level + 1);
623 623
            } else {
624 624
              _level->liftHighestActiveToTop();
625 625
            }
626 626
            if (_level->emptyLevel(level)) {
627 627
              _level->liftToTop(level);
628 628
            }
629 629
          } else {
630 630
            _level->deactivate(n);
631 631
          }
632 632

	
633 633
          n = _level->highestActive();
634 634
          level = _level->highestActiveLevel();
635 635
          --num;
636 636
        }
637 637

	
638 638
        num = _node_num * 20;
639 639
        while (num > 0 && n != INVALID) {
640
          Value excess = (*_excess)[n];
640
          Flow excess = (*_excess)[n];
641 641
          int new_level = _level->maxLevel();
642 642

	
643 643
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
644
            Value rem = (*_capacity)[e] - (*_flow)[e];
644
            Flow rem = (*_capacity)[e] - (*_flow)[e];
645 645
            if (!_tolerance.positive(rem)) continue;
646 646
            Node v = _graph.target(e);
647 647
            if ((*_level)[v] < level) {
648 648
              if (!_level->active(v) && v != _target) {
649 649
                _level->activate(v);
650 650
              }
651 651
              if (!_tolerance.less(rem, excess)) {
652 652
                _flow->set(e, (*_flow)[e] + excess);
653 653
                (*_excess)[v] += excess;
654 654
                excess = 0;
655 655
                goto no_more_push_2;
656 656
              } else {
657 657
                excess -= rem;
658 658
                (*_excess)[v] += rem;
659 659
                _flow->set(e, (*_capacity)[e]);
660 660
              }
661 661
            } else if (new_level > (*_level)[v]) {
662 662
              new_level = (*_level)[v];
663 663
            }
664 664
          }
665 665

	
666 666
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
667
            Value rem = (*_flow)[e];
667
            Flow rem = (*_flow)[e];
668 668
            if (!_tolerance.positive(rem)) continue;
669 669
            Node v = _graph.source(e);
670 670
            if ((*_level)[v] < level) {
671 671
              if (!_level->active(v) && v != _target) {
672 672
                _level->activate(v);
673 673
              }
674 674
              if (!_tolerance.less(rem, excess)) {
675 675
                _flow->set(e, (*_flow)[e] - excess);
676 676
                (*_excess)[v] += excess;
677 677
                excess = 0;
678 678
                goto no_more_push_2;
679 679
              } else {
680 680
                excess -= rem;
681 681
                (*_excess)[v] += rem;
682 682
                _flow->set(e, 0);
683 683
              }
684 684
            } else if (new_level > (*_level)[v]) {
685 685
              new_level = (*_level)[v];
686 686
            }
687 687
          }
688 688

	
689 689
        no_more_push_2:
690 690

	
691 691
          (*_excess)[n] = excess;
692 692

	
693 693
          if (excess != 0) {
694 694
            if (new_level + 1 < _level->maxLevel()) {
695 695
              _level->liftActiveOn(level, new_level + 1);
696 696
            } else {
697 697
              _level->liftActiveToTop(level);
698 698
            }
699 699
            if (_level->emptyLevel(level)) {
700 700
              _level->liftToTop(level);
701 701
            }
702 702
          } else {
703 703
            _level->deactivate(n);
704 704
          }
705 705

	
706 706
          while (level >= 0 && _level->activeFree(level)) {
707 707
            --level;
708 708
          }
709 709
          if (level == -1) {
710 710
            n = _level->highestActive();
711 711
            level = _level->highestActiveLevel();
712 712
          } else {
713 713
            n = _level->activeOn(level);
714 714
          }
715 715
          --num;
716 716
        }
717 717
      }
718 718
    }
719 719

	
720 720
    /// \brief Starts the second phase of the preflow algorithm.
721 721
    ///
722 722
    /// The preflow algorithm consists of two phases, this method runs
723 723
    /// the second phase. After calling one of the \ref init() functions
724 724
    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
725 725
    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
726 726
    /// value of a maximum flow, \ref minCut() returns a minimum cut
727 727
    /// \pre One of the \ref init() functions and \ref startFirstPhase()
728 728
    /// must be called before using this function.
729 729
    void startSecondPhase() {
730 730
      _phase = false;
731 731

	
732 732
      typename Digraph::template NodeMap<bool> reached(_graph);
733 733
      for (NodeIt n(_graph); n != INVALID; ++n) {
734 734
        reached[n] = (*_level)[n] < _level->maxLevel();
735 735
      }
736 736

	
737 737
      _level->initStart();
738 738
      _level->initAddItem(_source);
739 739

	
740 740
      std::vector<Node> queue;
741 741
      queue.push_back(_source);
742 742
      reached[_source] = true;
743 743

	
744 744
      while (!queue.empty()) {
745 745
        _level->initNewLevel();
746 746
        std::vector<Node> nqueue;
747 747
        for (int i = 0; i < int(queue.size()); ++i) {
748 748
          Node n = queue[i];
749 749
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
750 750
            Node v = _graph.target(e);
751 751
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
752 752
              reached[v] = true;
753 753
              _level->initAddItem(v);
754 754
              nqueue.push_back(v);
755 755
            }
756 756
          }
757 757
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
758 758
            Node u = _graph.source(e);
759 759
            if (!reached[u] &&
760 760
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
761 761
              reached[u] = true;
762 762
              _level->initAddItem(u);
763 763
              nqueue.push_back(u);
764 764
            }
765 765
          }
766 766
        }
767 767
        queue.swap(nqueue);
768 768
      }
769 769
      _level->initFinish();
770 770

	
771 771
      for (NodeIt n(_graph); n != INVALID; ++n) {
772 772
        if (!reached[n]) {
773 773
          _level->dirtyTopButOne(n);
774 774
        } else if ((*_excess)[n] > 0 && _target != n) {
775 775
          _level->activate(n);
776 776
        }
777 777
      }
778 778

	
779 779
      Node n;
780 780
      while ((n = _level->highestActive()) != INVALID) {
781
        Value excess = (*_excess)[n];
781
        Flow excess = (*_excess)[n];
782 782
        int level = _level->highestActiveLevel();
783 783
        int new_level = _level->maxLevel();
784 784

	
785 785
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
786
          Value rem = (*_capacity)[e] - (*_flow)[e];
786
          Flow rem = (*_capacity)[e] - (*_flow)[e];
787 787
          if (!_tolerance.positive(rem)) continue;
788 788
          Node v = _graph.target(e);
789 789
          if ((*_level)[v] < level) {
790 790
            if (!_level->active(v) && v != _source) {
791 791
              _level->activate(v);
792 792
            }
793 793
            if (!_tolerance.less(rem, excess)) {
794 794
              _flow->set(e, (*_flow)[e] + excess);
795 795
              (*_excess)[v] += excess;
796 796
              excess = 0;
797 797
              goto no_more_push;
798 798
            } else {
799 799
              excess -= rem;
800 800
              (*_excess)[v] += rem;
801 801
              _flow->set(e, (*_capacity)[e]);
802 802
            }
803 803
          } else if (new_level > (*_level)[v]) {
804 804
            new_level = (*_level)[v];
805 805
          }
806 806
        }
807 807

	
808 808
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
809
          Value rem = (*_flow)[e];
809
          Flow rem = (*_flow)[e];
810 810
          if (!_tolerance.positive(rem)) continue;
811 811
          Node v = _graph.source(e);
812 812
          if ((*_level)[v] < level) {
813 813
            if (!_level->active(v) && v != _source) {
814 814
              _level->activate(v);
815 815
            }
816 816
            if (!_tolerance.less(rem, excess)) {
817 817
              _flow->set(e, (*_flow)[e] - excess);
818 818
              (*_excess)[v] += excess;
819 819
              excess = 0;
820 820
              goto no_more_push;
821 821
            } else {
822 822
              excess -= rem;
823 823
              (*_excess)[v] += rem;
824 824
              _flow->set(e, 0);
825 825
            }
826 826
          } else if (new_level > (*_level)[v]) {
827 827
            new_level = (*_level)[v];
828 828
          }
829 829
        }
830 830

	
831 831
      no_more_push:
832 832

	
833 833
        (*_excess)[n] = excess;
834 834

	
835 835
        if (excess != 0) {
836 836
          if (new_level + 1 < _level->maxLevel()) {
837 837
            _level->liftHighestActive(new_level + 1);
838 838
          } else {
839 839
            // Calculation error
840 840
            _level->liftHighestActiveToTop();
841 841
          }
842 842
          if (_level->emptyLevel(level)) {
843 843
            // Calculation error
844 844
            _level->liftToTop(level);
845 845
          }
846 846
        } else {
847 847
          _level->deactivate(n);
848 848
        }
849 849

	
850 850
      }
851 851
    }
852 852

	
853 853
    /// \brief Runs the preflow algorithm.
854 854
    ///
855 855
    /// Runs the preflow algorithm.
856 856
    /// \note pf.run() is just a shortcut of the following code.
857 857
    /// \code
858 858
    ///   pf.init();
859 859
    ///   pf.startFirstPhase();
860 860
    ///   pf.startSecondPhase();
861 861
    /// \endcode
862 862
    void run() {
863 863
      init();
864 864
      startFirstPhase();
865 865
      startSecondPhase();
866 866
    }
867 867

	
868 868
    /// \brief Runs the preflow algorithm to compute the minimum cut.
869 869
    ///
870 870
    /// Runs the preflow algorithm to compute the minimum cut.
871 871
    /// \note pf.runMinCut() is just a shortcut of the following code.
872 872
    /// \code
873 873
    ///   pf.init();
874 874
    ///   pf.startFirstPhase();
875 875
    /// \endcode
876 876
    void runMinCut() {
877 877
      init();
878 878
      startFirstPhase();
879 879
    }
880 880

	
881 881
    /// @}
882 882

	
883 883
    /// \name Query Functions
884 884
    /// The results of the preflow algorithm can be obtained using these
885 885
    /// functions.\n
886 886
    /// Either one of the \ref run() "run*()" functions or one of the
887 887
    /// \ref startFirstPhase() "start*()" functions should be called
888 888
    /// before using them.
889 889

	
890 890
    ///@{
891 891

	
892 892
    /// \brief Returns the value of the maximum flow.
893 893
    ///
894 894
    /// Returns the value of the maximum flow by returning the excess
895 895
    /// of the target node. This value equals to the value of
896 896
    /// the maximum flow already after the first phase of the algorithm.
897 897
    ///
898 898
    /// \pre Either \ref run() or \ref init() must be called before
899 899
    /// using this function.
900
    Value flowValue() const {
900
    Flow flowValue() const {
901 901
      return (*_excess)[_target];
902 902
    }
903 903

	
904 904
    /// \brief Returns the flow on the given arc.
905 905
    ///
906 906
    /// Returns the flow on the given arc. This method can
907 907
    /// be called after the second phase of the algorithm.
908 908
    ///
909 909
    /// \pre Either \ref run() or \ref init() must be called before
910 910
    /// using this function.
911
    Value flow(const Arc& arc) const {
911
    Flow flow(const Arc& arc) const {
912 912
      return (*_flow)[arc];
913 913
    }
914 914

	
915 915
    /// \brief Returns a const reference to the flow map.
916 916
    ///
917 917
    /// Returns a const reference to the arc map storing the found flow.
918 918
    /// This method can be called after the second phase of the algorithm.
919 919
    ///
920 920
    /// \pre Either \ref run() or \ref init() must be called before
921 921
    /// using this function.
922 922
    const FlowMap& flowMap() const {
923 923
      return *_flow;
924 924
    }
925 925

	
926 926
    /// \brief Returns \c true when the node is on the source side of the
927 927
    /// minimum cut.
928 928
    ///
929 929
    /// Returns true when the node is on the source side of the found
930 930
    /// minimum cut. This method can be called both after running \ref
931 931
    /// startFirstPhase() and \ref startSecondPhase().
932 932
    ///
933 933
    /// \pre Either \ref run() or \ref init() must be called before
934 934
    /// using this function.
935 935
    bool minCut(const Node& node) const {
936 936
      return ((*_level)[node] == _level->maxLevel()) == _phase;
937 937
    }
938 938

	
939 939
    /// \brief Gives back a minimum value cut.
940 940
    ///
941 941
    /// Sets \c cutMap to the characteristic vector of a minimum value
942 942
    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
943 943
    /// node map with \c bool (or convertible) value type.
944 944
    ///
945 945
    /// This method can be called both after running \ref startFirstPhase()
946 946
    /// and \ref startSecondPhase(). The result after the second phase
947 947
    /// could be slightly different if inexact computation is used.
948 948
    ///
949 949
    /// \note This function calls \ref minCut() for each node, so it runs in
950 950
    /// O(n) time.
951 951
    ///
952 952
    /// \pre Either \ref run() or \ref init() must be called before
953 953
    /// using this function.
954 954
    template <typename CutMap>
955 955
    void minCutMap(CutMap& cutMap) const {
956 956
      for (NodeIt n(_graph); n != INVALID; ++n) {
957 957
        cutMap.set(n, minCut(n));
958 958
      }
959 959
    }
960 960

	
961 961
    /// @}
962 962
  };
963 963
}
964 964

	
965 965
#endif
Ignore white space 6 line context
1 1
INCLUDE_DIRECTORIES(
2 2
  ${PROJECT_SOURCE_DIR}
3 3
  ${PROJECT_BINARY_DIR}
4 4
)
5 5

	
6 6
IF(HAVE_GLPK)
7 7
  INCLUDE_DIRECTORIES(${GLPK_INCLUDE_DIR})
8 8
ENDIF(HAVE_GLPK)
9 9

	
10 10
LINK_DIRECTORIES(${PROJECT_BINARY_DIR}/lemon)
11 11

	
12 12
SET(TESTS
13 13
  adaptors_test
14 14
  bfs_test
15 15
  circulation_test
16 16
  counter_test
17 17
  dfs_test
18 18
  digraph_test
19 19
  dijkstra_test
20 20
  dim_test
21 21
  edge_set_test
22 22
  error_test
23 23
  euler_test
24 24
  gomory_hu_test
25 25
  graph_copy_test
26 26
  graph_test
27 27
  graph_utils_test
28 28
  hao_orlin_test
29 29
  heap_test
30 30
  kruskal_test
31 31
  maps_test
32 32
  matching_test
33 33
  min_cost_arborescence_test
34
  min_cost_flow_test
34 35
  path_test
35 36
  preflow_test
36 37
  radix_sort_test
37 38
  random_test
38 39
  suurballe_test
39 40
  time_measure_test
40 41
  unionfind_test)
41 42

	
42 43
IF(HAVE_LP)
43 44
  ADD_EXECUTABLE(lp_test lp_test.cc)
44 45
  IF(HAVE_GLPK)
45 46
    TARGET_LINK_LIBRARIES(lp_test lemon ${GLPK_LIBRARIES})
46 47
  ENDIF(HAVE_GLPK)
47 48
  ADD_TEST(lp_test lp_test)
48 49

	
49 50
  IF(WIN32 AND HAVE_GLPK)
50 51
    GET_TARGET_PROPERTY(TARGET_LOC lp_test LOCATION)
51 52
    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
52 53
    ADD_CUSTOM_COMMAND(TARGET lp_test POST_BUILD
53 54
      COMMAND cmake -E copy ${GLPK_BIN_DIR}/glpk.dll ${TARGET_PATH}
54 55
      COMMAND cmake -E copy ${GLPK_BIN_DIR}/libltdl3.dll ${TARGET_PATH}
55 56
      COMMAND cmake -E copy ${GLPK_BIN_DIR}/zlib1.dll ${TARGET_PATH}
56 57
    )
57 58
  ENDIF(WIN32 AND HAVE_GLPK)
58 59
ENDIF(HAVE_LP)
59 60

	
60 61
IF(HAVE_MIP)
61 62
  ADD_EXECUTABLE(mip_test mip_test.cc)
62 63
  IF(HAVE_GLPK)
63 64
    TARGET_LINK_LIBRARIES(mip_test lemon ${GLPK_LIBRARIES})
64 65
  ENDIF(HAVE_GLPK)
65 66
  ADD_TEST(mip_test mip_test)
66 67

	
67 68
  IF(WIN32 AND HAVE_GLPK)
68 69
    GET_TARGET_PROPERTY(TARGET_LOC mip_test LOCATION)
69 70
    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
70 71
    ADD_CUSTOM_COMMAND(TARGET mip_test POST_BUILD
71 72
      COMMAND cmake -E copy ${GLPK_BIN_DIR}/glpk.dll ${TARGET_PATH}
72 73
      COMMAND cmake -E copy ${GLPK_BIN_DIR}/libltdl3.dll ${TARGET_PATH}
73 74
      COMMAND cmake -E copy ${GLPK_BIN_DIR}/zlib1.dll ${TARGET_PATH}
74 75
    )
75 76
  ENDIF(WIN32 AND HAVE_GLPK)
76 77
ENDIF(HAVE_MIP)
77 78

	
78 79
FOREACH(TEST_NAME ${TESTS})
79 80
  ADD_EXECUTABLE(${TEST_NAME} ${TEST_NAME}.cc)
80 81
  TARGET_LINK_LIBRARIES(${TEST_NAME} lemon)
81 82
  ADD_TEST(${TEST_NAME} ${TEST_NAME})
82 83
ENDFOREACH(TEST_NAME)
Ignore white space 6 line context
1 1
EXTRA_DIST += \
2 2
	test/CMakeLists.txt
3 3

	
4 4
noinst_HEADERS += \
5 5
	test/graph_test.h \
6 6
	test/test_tools.h
7 7

	
8 8
check_PROGRAMS += \
9 9
	test/adaptors_test \
10 10
	test/bfs_test \
11 11
	test/circulation_test \
12 12
	test/counter_test \
13 13
	test/dfs_test \
14 14
	test/digraph_test \
15 15
	test/dijkstra_test \
16 16
	test/dim_test \
17 17
	test/edge_set_test \
18 18
	test/error_test \
19 19
	test/euler_test \
20 20
	test/gomory_hu_test \
21 21
	test/graph_copy_test \
22 22
	test/graph_test \
23 23
	test/graph_utils_test \
24 24
	test/hao_orlin_test \
25 25
	test/heap_test \
26 26
	test/kruskal_test \
27 27
	test/maps_test \
28 28
	test/matching_test \
29 29
	test/min_cost_arborescence_test \
30
	test/min_cost_flow_test \
30 31
	test/path_test \
31 32
	test/preflow_test \
32 33
	test/radix_sort_test \
33 34
	test/random_test \
34 35
	test/suurballe_test \
35 36
	test/test_tools_fail \
36 37
	test/test_tools_pass \
37 38
	test/time_measure_test \
38 39
	test/unionfind_test
39 40

	
40 41
test_test_tools_pass_DEPENDENCIES = demo
41 42

	
42 43
if HAVE_LP
43 44
check_PROGRAMS += test/lp_test
44 45
endif HAVE_LP
45 46
if HAVE_MIP
46 47
check_PROGRAMS += test/mip_test
47 48
endif HAVE_MIP
48 49

	
49 50
TESTS += $(check_PROGRAMS)
50 51
XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
51 52

	
52 53
test_adaptors_test_SOURCES = test/adaptors_test.cc
53 54
test_bfs_test_SOURCES = test/bfs_test.cc
54 55
test_circulation_test_SOURCES = test/circulation_test.cc
55 56
test_counter_test_SOURCES = test/counter_test.cc
56 57
test_dfs_test_SOURCES = test/dfs_test.cc
57 58
test_digraph_test_SOURCES = test/digraph_test.cc
58 59
test_dijkstra_test_SOURCES = test/dijkstra_test.cc
59 60
test_dim_test_SOURCES = test/dim_test.cc
60 61
test_edge_set_test_SOURCES = test/edge_set_test.cc
61 62
test_error_test_SOURCES = test/error_test.cc
62 63
test_euler_test_SOURCES = test/euler_test.cc
63 64
test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
64 65
test_graph_copy_test_SOURCES = test/graph_copy_test.cc
65 66
test_graph_test_SOURCES = test/graph_test.cc
66 67
test_graph_utils_test_SOURCES = test/graph_utils_test.cc
67 68
test_heap_test_SOURCES = test/heap_test.cc
68 69
test_kruskal_test_SOURCES = test/kruskal_test.cc
69 70
test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
70 71
test_lp_test_SOURCES = test/lp_test.cc
71 72
test_maps_test_SOURCES = test/maps_test.cc
72 73
test_mip_test_SOURCES = test/mip_test.cc
73 74
test_matching_test_SOURCES = test/matching_test.cc
74 75
test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
76
test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
75 77
test_path_test_SOURCES = test/path_test.cc
76 78
test_preflow_test_SOURCES = test/preflow_test.cc
77 79
test_radix_sort_test_SOURCES = test/radix_sort_test.cc
78 80
test_suurballe_test_SOURCES = test/suurballe_test.cc
79 81
test_random_test_SOURCES = test/random_test.cc
80 82
test_test_tools_fail_SOURCES = test/test_tools_fail.cc
81 83
test_test_tools_pass_SOURCES = test/test_tools_pass.cc
82 84
test_time_measure_test_SOURCES = test/time_measure_test.cc
83 85
test_unionfind_test_SOURCES = test/unionfind_test.cc
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#include <iostream>
20 20

	
21 21
#include "test_tools.h"
22 22
#include <lemon/list_graph.h>
23 23
#include <lemon/circulation.h>
24 24
#include <lemon/lgf_reader.h>
25 25
#include <lemon/concepts/digraph.h>
26 26
#include <lemon/concepts/maps.h>
27 27

	
28 28
using namespace lemon;
29 29

	
30 30
char test_lgf[] =
31 31
  "@nodes\n"
32 32
  "label\n"
33 33
  "0\n"
34 34
  "1\n"
35 35
  "2\n"
36 36
  "3\n"
37 37
  "4\n"
38 38
  "5\n"
39 39
  "@arcs\n"
40 40
  "     lcap  ucap\n"
41 41
  "0 1  2  10\n"
42 42
  "0 2  2  6\n"
43 43
  "1 3  4  7\n"
44 44
  "1 4  0  5\n"
45 45
  "2 4  1  3\n"
46 46
  "3 5  3  8\n"
47 47
  "4 5  3  7\n"
48 48
  "@attributes\n"
49 49
  "source 0\n"
50 50
  "sink   5\n";
51 51

	
52 52
void checkCirculationCompile()
53 53
{
54 54
  typedef int VType;
55 55
  typedef concepts::Digraph Digraph;
56 56

	
57 57
  typedef Digraph::Node Node;
58 58
  typedef Digraph::Arc Arc;
59 59
  typedef concepts::ReadMap<Arc,VType> CapMap;
60
  typedef concepts::ReadMap<Node,VType> DeltaMap;
60
  typedef concepts::ReadMap<Node,VType> SupplyMap;
61 61
  typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
62 62
  typedef concepts::WriteMap<Node,bool> BarrierMap;
63 63

	
64 64
  typedef Elevator<Digraph, Digraph::Node> Elev;
65 65
  typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev;
66 66

	
67 67
  Digraph g;
68 68
  Node n;
69 69
  Arc a;
70 70
  CapMap lcap, ucap;
71
  DeltaMap delta;
71
  SupplyMap supply;
72 72
  FlowMap flow;
73 73
  BarrierMap bar;
74 74
  VType v;
75 75
  bool b;
76 76

	
77
  typedef Circulation<Digraph, CapMap, CapMap, DeltaMap>
77
  typedef Circulation<Digraph, CapMap, CapMap, SupplyMap>
78 78
            ::SetFlowMap<FlowMap>
79 79
            ::SetElevator<Elev>
80 80
            ::SetStandardElevator<LinkedElev>
81 81
            ::Create CirculationType;
82
  CirculationType circ_test(g, lcap, ucap, delta);
82
  CirculationType circ_test(g, lcap, ucap, supply);
83 83
  const CirculationType& const_circ_test = circ_test;
84 84
   
85 85
  circ_test
86
    .lowerCapMap(lcap)
87
    .upperCapMap(ucap)
88
    .deltaMap(delta)
86
    .lowerMap(lcap)
87
    .upperMap(ucap)
88
    .supplyMap(supply)
89 89
    .flowMap(flow);
90 90

	
91 91
  circ_test.init();
92 92
  circ_test.greedyInit();
93 93
  circ_test.start();
94 94
  circ_test.run();
95 95

	
96 96
  v = const_circ_test.flow(a);
97 97
  const FlowMap& fm = const_circ_test.flowMap();
98 98
  b = const_circ_test.barrier(n);
99 99
  const_circ_test.barrierMap(bar);
100 100
  
101 101
  ignore_unused_variable_warning(fm);
102 102
}
103 103

	
104 104
template <class G, class LM, class UM, class DM>
105 105
void checkCirculation(const G& g, const LM& lm, const UM& um,
106 106
                      const DM& dm, bool find)
107 107
{
108 108
  Circulation<G, LM, UM, DM> circ(g, lm, um, dm);
109 109
  bool ret = circ.run();
110 110
  if (find) {
111 111
    check(ret, "A feasible solution should have been found.");
112 112
    check(circ.checkFlow(), "The found flow is corrupt.");
113 113
    check(!circ.checkBarrier(), "A barrier should not have been found.");
114 114
  } else {
115 115
    check(!ret, "A feasible solution should not have been found.");
116 116
    check(circ.checkBarrier(), "The found barrier is corrupt.");
117 117
  }
118 118
}
119 119

	
120 120
int main (int, char*[])
121 121
{
122 122
  typedef ListDigraph Digraph;
123 123
  DIGRAPH_TYPEDEFS(Digraph);
124 124

	
125 125
  Digraph g;
126 126
  IntArcMap lo(g), up(g);
127 127
  IntNodeMap delta(g, 0);
128 128
  Node s, t;
129 129

	
130 130
  std::istringstream input(test_lgf);
131 131
  DigraphReader<Digraph>(g,input).
132 132
    arcMap("lcap", lo).
133 133
    arcMap("ucap", up).
134 134
    node("source",s).
135 135
    node("sink",t).
136 136
    run();
137 137

	
138 138
  delta[s] = 7; delta[t] = -7;
139 139
  checkCirculation(g, lo, up, delta, true);
140 140

	
141 141
  delta[s] = 13; delta[t] = -13;
142 142
  checkCirculation(g, lo, up, delta, true);
143 143

	
144 144
  delta[s] = 6; delta[t] = -6;
145 145
  checkCirculation(g, lo, up, delta, false);
146 146

	
147 147
  delta[s] = 14; delta[t] = -14;
148 148
  checkCirculation(g, lo, up, delta, false);
149 149

	
150 150
  delta[s] = 7; delta[t] = -13;
151 151
  checkCirculation(g, lo, up, delta, true);
152 152

	
153 153
  delta[s] = 5; delta[t] = -15;
154 154
  checkCirculation(g, lo, up, delta, true);
155 155

	
156 156
  delta[s] = 10; delta[t] = -11;
157 157
  checkCirculation(g, lo, up, delta, true);
158 158

	
159 159
  delta[s] = 11; delta[t] = -10;
160 160
  checkCirculation(g, lo, up, delta, false);
161 161

	
162 162
  return 0;
163 163
}
Ignore white space 3072 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
///\ingroup tools
20 20
///\file
21 21
///\brief DIMACS problem solver.
22 22
///
23 23
/// This program solves various problems given in DIMACS format.
24 24
///
25 25
/// See
26 26
/// \code
27 27
///   dimacs-solver --help
28 28
/// \endcode
29 29
/// for more info on usage.
30 30

	
31 31
#include <iostream>
32 32
#include <fstream>
33 33
#include <cstring>
34 34

	
35 35
#include <lemon/smart_graph.h>
36 36
#include <lemon/dimacs.h>
37 37
#include <lemon/lgf_writer.h>
38 38
#include <lemon/time_measure.h>
39 39

	
40 40
#include <lemon/arg_parser.h>
41 41
#include <lemon/error.h>
42 42

	
43 43
#include <lemon/dijkstra.h>
44 44
#include <lemon/preflow.h>
45 45
#include <lemon/matching.h>
46
#include <lemon/network_simplex.h>
46 47

	
47 48
using namespace lemon;
48 49
typedef SmartDigraph Digraph;
49 50
DIGRAPH_TYPEDEFS(Digraph);
50 51
typedef SmartGraph Graph;
51 52

	
52 53
template<class Value>
53 54
void solve_sp(ArgParser &ap, std::istream &is, std::ostream &,
54 55
              DimacsDescriptor &desc)
55 56
{
56 57
  bool report = !ap.given("q");
57 58
  Digraph g;
58 59
  Node s;
59 60
  Digraph::ArcMap<Value> len(g);
60 61
  Timer t;
61 62
  t.restart();
62 63
  readDimacsSp(is, g, len, s, desc);
63 64
  if(report) std::cerr << "Read the file: " << t << '\n';
64 65
  t.restart();
65 66
  Dijkstra<Digraph, Digraph::ArcMap<Value> > dij(g,len);
66 67
  if(report) std::cerr << "Setup Dijkstra class: " << t << '\n';
67 68
  t.restart();
68 69
  dij.run(s);
69 70
  if(report) std::cerr << "Run Dijkstra: " << t << '\n';
70 71
}
71 72

	
72 73
template<class Value>
73 74
void solve_max(ArgParser &ap, std::istream &is, std::ostream &,
74 75
               Value infty, DimacsDescriptor &desc)
75 76
{
76 77
  bool report = !ap.given("q");
77 78
  Digraph g;
78 79
  Node s,t;
79 80
  Digraph::ArcMap<Value> cap(g);
80 81
  Timer ti;
81 82
  ti.restart();
82 83
  readDimacsMax(is, g, cap, s, t, infty, desc);
83 84
  if(report) std::cerr << "Read the file: " << ti << '\n';
84 85
  ti.restart();
85 86
  Preflow<Digraph, Digraph::ArcMap<Value> > pre(g,cap,s,t);
86 87
  if(report) std::cerr << "Setup Preflow class: " << ti << '\n';
87 88
  ti.restart();
88 89
  pre.run();
89 90
  if(report) std::cerr << "Run Preflow: " << ti << '\n';
90 91
  if(report) std::cerr << "\nMax flow value: " << pre.flowValue() << '\n';  
91 92
}
92 93

	
94
template<class Value>
95
void solve_min(ArgParser &ap, std::istream &is, std::ostream &,
96
               DimacsDescriptor &desc)
97
{
98
  bool report = !ap.given("q");
99
  Digraph g;
100
  Digraph::ArcMap<Value> lower(g), cap(g), cost(g);
101
  Digraph::NodeMap<Value> sup(g);
102
  Timer ti;
103
  ti.restart();
104
  readDimacsMin(is, g, lower, cap, cost, sup, 0, desc);
105
  if (report) std::cerr << "Read the file: " << ti << '\n';
106
  ti.restart();
107
  NetworkSimplex<Digraph, Value> ns(g);
108
  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
109
  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
110
  ti.restart();
111
  ns.run();
112
  if (report) std::cerr << "Run NetworkSimplex: " << ti << '\n';
113
  if (report) std::cerr << "\nMin flow cost: " << ns.totalCost() << '\n';
114
}
115

	
93 116
void solve_mat(ArgParser &ap, std::istream &is, std::ostream &,
94 117
              DimacsDescriptor &desc)
95 118
{
96 119
  bool report = !ap.given("q");
97 120
  Graph g;
98 121
  Timer ti;
99 122
  ti.restart();
100 123
  readDimacsMat(is, g, desc);
101 124
  if(report) std::cerr << "Read the file: " << ti << '\n';
102 125
  ti.restart();
103 126
  MaxMatching<Graph> mat(g);
104 127
  if(report) std::cerr << "Setup MaxMatching class: " << ti << '\n';
105 128
  ti.restart();
106 129
  mat.run();
107 130
  if(report) std::cerr << "Run MaxMatching: " << ti << '\n';
108 131
  if(report) std::cerr << "\nCardinality of max matching: "
109 132
                       << mat.matchingSize() << '\n';  
110 133
}
111 134

	
112 135

	
113 136
template<class Value>
114 137
void solve(ArgParser &ap, std::istream &is, std::ostream &os,
115 138
           DimacsDescriptor &desc)
116 139
{
117 140
  std::stringstream iss(static_cast<std::string>(ap["infcap"]));
118 141
  Value infty;
119 142
  iss >> infty;
120 143
  if(iss.fail())
121 144
    {
122 145
      std::cerr << "Cannot interpret '"
123 146
                << static_cast<std::string>(ap["infcap"]) << "' as infinite"
124 147
                << std::endl;
125 148
      exit(1);
126 149
    }
127 150
  
128 151
  switch(desc.type)
129 152
    {
130 153
    case DimacsDescriptor::MIN:
131
      std::cerr <<
132
        "\n\n Sorry, the min. cost flow solver is not yet available.\n";
154
      solve_min<Value>(ap,is,os,desc);
133 155
      break;
134 156
    case DimacsDescriptor::MAX:
135 157
      solve_max<Value>(ap,is,os,infty,desc);
136 158
      break;
137 159
    case DimacsDescriptor::SP:
138 160
      solve_sp<Value>(ap,is,os,desc);
139 161
      break;
140 162
    case DimacsDescriptor::MAT:
141 163
      solve_mat(ap,is,os,desc);
142 164
      break;
143 165
    default:
144 166
      break;
145 167
    }
146 168
}
147 169

	
148 170
int main(int argc, const char *argv[]) {
149 171
  typedef SmartDigraph Digraph;
150 172

	
151 173
  typedef Digraph::Arc Arc;
152 174

	
153 175
  std::string inputName;
154 176
  std::string outputName;
155 177

	
156 178
  ArgParser ap(argc, argv);
157 179
  ap.other("[INFILE [OUTFILE]]",
158 180
           "If either the INFILE or OUTFILE file is missing the standard\n"
159 181
           "     input/output will be used instead.")
160 182
    .boolOption("q", "Do not print any report")
161 183
    .boolOption("int","Use 'int' for capacities, costs etc. (default)")
162 184
    .optionGroup("datatype","int")
163 185
#ifdef HAVE_LONG_LONG
164 186
    .boolOption("long","Use 'long long' for capacities, costs etc.")
165 187
    .optionGroup("datatype","long")
166 188
#endif
167 189
    .boolOption("double","Use 'double' for capacities, costs etc.")
168 190
    .optionGroup("datatype","double")
169 191
    .boolOption("ldouble","Use 'long double' for capacities, costs etc.")
170 192
    .optionGroup("datatype","ldouble")
171 193
    .onlyOneGroup("datatype")
172 194
    .stringOption("infcap","Value used for 'very high' capacities","0")
173 195
    .run();
174 196

	
175 197
  std::ifstream input;
176 198
  std::ofstream output;
177 199

	
178 200
  switch(ap.files().size())
179 201
    {
180 202
    case 2:
181 203
      output.open(ap.files()[1].c_str());
182 204
      if (!output) {
183 205
        throw IoError("Cannot open the file for writing", ap.files()[1]);
184 206
      }
185 207
    case 1:
186 208
      input.open(ap.files()[0].c_str());
187 209
      if (!input) {
188 210
        throw IoError("File cannot be found", ap.files()[0]);
189 211
      }
190 212
    case 0:
191 213
      break;
192 214
    default:
193 215
      std::cerr << ap.commandName() << ": too many arguments\n";
194 216
      return 1;
195 217
    }
196 218
  std::istream& is = (ap.files().size()<1 ? std::cin : input);
197 219
  std::ostream& os = (ap.files().size()<2 ? std::cout : output);
198 220

	
199 221
  DimacsDescriptor desc = dimacsType(is);
200 222
  
201 223
  if(!ap.given("q"))
202 224
    {
203 225
      std::cout << "Problem type: ";
204 226
      switch(desc.type)
205 227
        {
206 228
        case DimacsDescriptor::MIN:
207 229
          std::cout << "min";
208 230
          break;
209 231
        case DimacsDescriptor::MAX:
210 232
          std::cout << "max";
211 233
          break;
212 234
        case DimacsDescriptor::SP:
213 235
          std::cout << "sp";
214 236
        case DimacsDescriptor::MAT:
215 237
          std::cout << "mat";
216 238
          break;
217 239
        default:
218 240
          exit(1);
219 241
          break;
220 242
        }
221 243
      std::cout << "\nNum of nodes: " << desc.nodeNum;
222 244
      std::cout << "\nNum of arcs:  " << desc.edgeNum;
223 245
      std::cout << "\n\n";
224 246
    }
225 247
    
226 248
  if(ap.given("double"))
227 249
    solve<double>(ap,is,os,desc);
228 250
  else if(ap.given("ldouble"))
229 251
    solve<long double>(ap,is,os,desc);
230 252
#ifdef HAVE_LONG_LONG
231 253
  else if(ap.given("long"))
232 254
    solve<long long>(ap,is,os,desc);
233 255
#endif
234 256
  else solve<int>(ap,is,os,desc);
235 257

	
236 258
  return 0;
237 259
}
0 comments (0 inline)