lemon/network_simplex.h
author Peter Kovacs <kpeter@inf.elte.hu>
Tue, 24 Feb 2009 09:46:02 +0100
changeset 601 e8349c6f12ca
child 603 425cc8328c0e
permissions -rw-r--r--
Port NetworkSimplex from SVN -r3520 (#234)
kpeter@601
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
kpeter@601
     2
 *
kpeter@601
     3
 * This file is a part of LEMON, a generic C++ optimization library.
kpeter@601
     4
 *
kpeter@601
     5
 * Copyright (C) 2003-2009
kpeter@601
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
kpeter@601
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
kpeter@601
     8
 *
kpeter@601
     9
 * Permission to use, modify and distribute this software is granted
kpeter@601
    10
 * provided that this copyright notice appears in all copies. For
kpeter@601
    11
 * precise terms see the accompanying LICENSE file.
kpeter@601
    12
 *
kpeter@601
    13
 * This software is provided "AS IS" with no warranty of any kind,
kpeter@601
    14
 * express or implied, and with no claim as to its suitability for any
kpeter@601
    15
 * purpose.
kpeter@601
    16
 *
kpeter@601
    17
 */
kpeter@601
    18
kpeter@601
    19
#ifndef LEMON_NETWORK_SIMPLEX_H
kpeter@601
    20
#define LEMON_NETWORK_SIMPLEX_H
kpeter@601
    21
kpeter@601
    22
/// \ingroup min_cost_flow
kpeter@601
    23
///
kpeter@601
    24
/// \file
kpeter@601
    25
/// \brief Network simplex algorithm for finding a minimum cost flow.
kpeter@601
    26
kpeter@601
    27
#include <vector>
kpeter@601
    28
#include <limits>
kpeter@601
    29
#include <algorithm>
kpeter@601
    30
kpeter@601
    31
#include <lemon/math.h>
kpeter@601
    32
kpeter@601
    33
namespace lemon {
kpeter@601
    34
kpeter@601
    35
  /// \addtogroup min_cost_flow
kpeter@601
    36
  /// @{
kpeter@601
    37
kpeter@601
    38
  /// \brief Implementation of the primal network simplex algorithm
kpeter@601
    39
  /// for finding a \ref min_cost_flow "minimum cost flow".
kpeter@601
    40
  ///
kpeter@601
    41
  /// \ref NetworkSimplex implements the primal network simplex algorithm
kpeter@601
    42
  /// for finding a \ref min_cost_flow "minimum cost flow".
kpeter@601
    43
  ///
kpeter@601
    44
  /// \tparam Digraph The digraph type the algorithm runs on.
kpeter@601
    45
  /// \tparam LowerMap The type of the lower bound map.
kpeter@601
    46
  /// \tparam CapacityMap The type of the capacity (upper bound) map.
kpeter@601
    47
  /// \tparam CostMap The type of the cost (length) map.
kpeter@601
    48
  /// \tparam SupplyMap The type of the supply map.
kpeter@601
    49
  ///
kpeter@601
    50
  /// \warning
kpeter@601
    51
  /// - Arc capacities and costs should be \e non-negative \e integers.
kpeter@601
    52
  /// - Supply values should be \e signed \e integers.
kpeter@601
    53
  /// - The value types of the maps should be convertible to each other.
kpeter@601
    54
  /// - \c CostMap::Value must be signed type.
kpeter@601
    55
  ///
kpeter@601
    56
  /// \note \ref NetworkSimplex provides five different pivot rule
kpeter@601
    57
  /// implementations that significantly affect the efficiency of the
kpeter@601
    58
  /// algorithm.
kpeter@601
    59
  /// By default "Block Search" pivot rule is used, which proved to be
kpeter@601
    60
  /// by far the most efficient according to our benchmark tests.
kpeter@601
    61
  /// However another pivot rule can be selected using \ref run()
kpeter@601
    62
  /// function with the proper parameter.
kpeter@601
    63
#ifdef DOXYGEN
kpeter@601
    64
  template < typename Digraph,
kpeter@601
    65
             typename LowerMap,
kpeter@601
    66
             typename CapacityMap,
kpeter@601
    67
             typename CostMap,
kpeter@601
    68
             typename SupplyMap >
kpeter@601
    69
kpeter@601
    70
#else
kpeter@601
    71
  template < typename Digraph,
kpeter@601
    72
             typename LowerMap = typename Digraph::template ArcMap<int>,
kpeter@601
    73
             typename CapacityMap = typename Digraph::template ArcMap<int>,
kpeter@601
    74
             typename CostMap = typename Digraph::template ArcMap<int>,
kpeter@601
    75
             typename SupplyMap = typename Digraph::template NodeMap<int> >
kpeter@601
    76
#endif
kpeter@601
    77
  class NetworkSimplex
kpeter@601
    78
  {
kpeter@601
    79
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
kpeter@601
    80
kpeter@601
    81
    typedef typename CapacityMap::Value Capacity;
kpeter@601
    82
    typedef typename CostMap::Value Cost;
kpeter@601
    83
    typedef typename SupplyMap::Value Supply;
kpeter@601
    84
kpeter@601
    85
    typedef std::vector<Arc> ArcVector;
kpeter@601
    86
    typedef std::vector<Node> NodeVector;
kpeter@601
    87
    typedef std::vector<int> IntVector;
kpeter@601
    88
    typedef std::vector<bool> BoolVector;
kpeter@601
    89
    typedef std::vector<Capacity> CapacityVector;
kpeter@601
    90
    typedef std::vector<Cost> CostVector;
kpeter@601
    91
    typedef std::vector<Supply> SupplyVector;
kpeter@601
    92
kpeter@601
    93
  public:
kpeter@601
    94
kpeter@601
    95
    /// The type of the flow map
kpeter@601
    96
    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
kpeter@601
    97
    /// The type of the potential map
kpeter@601
    98
    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
kpeter@601
    99
kpeter@601
   100
  public:
kpeter@601
   101
kpeter@601
   102
    /// Enum type for selecting the pivot rule used by \ref run()
kpeter@601
   103
    enum PivotRuleEnum {
kpeter@601
   104
      FIRST_ELIGIBLE_PIVOT,
kpeter@601
   105
      BEST_ELIGIBLE_PIVOT,
kpeter@601
   106
      BLOCK_SEARCH_PIVOT,
kpeter@601
   107
      CANDIDATE_LIST_PIVOT,
kpeter@601
   108
      ALTERING_LIST_PIVOT
kpeter@601
   109
    };
kpeter@601
   110
kpeter@601
   111
  private:
kpeter@601
   112
kpeter@601
   113
    // State constants for arcs
kpeter@601
   114
    enum ArcStateEnum {
kpeter@601
   115
      STATE_UPPER = -1,
kpeter@601
   116
      STATE_TREE  =  0,
kpeter@601
   117
      STATE_LOWER =  1
kpeter@601
   118
    };
kpeter@601
   119
kpeter@601
   120
  private:
kpeter@601
   121
kpeter@601
   122
    // References for the original data
kpeter@601
   123
    const Digraph &_orig_graph;
kpeter@601
   124
    const LowerMap *_orig_lower;
kpeter@601
   125
    const CapacityMap &_orig_cap;
kpeter@601
   126
    const CostMap &_orig_cost;
kpeter@601
   127
    const SupplyMap *_orig_supply;
kpeter@601
   128
    Node _orig_source;
kpeter@601
   129
    Node _orig_target;
kpeter@601
   130
    Capacity _orig_flow_value;
kpeter@601
   131
kpeter@601
   132
    // Result maps
kpeter@601
   133
    FlowMap *_flow_result;
kpeter@601
   134
    PotentialMap *_potential_result;
kpeter@601
   135
    bool _local_flow;
kpeter@601
   136
    bool _local_potential;
kpeter@601
   137
kpeter@601
   138
    // Data structures for storing the graph
kpeter@601
   139
    ArcVector _arc;
kpeter@601
   140
    NodeVector _node;
kpeter@601
   141
    IntNodeMap _node_id;
kpeter@601
   142
    IntVector _source;
kpeter@601
   143
    IntVector _target;
kpeter@601
   144
kpeter@601
   145
    // The number of nodes and arcs in the original graph
kpeter@601
   146
    int _node_num;
kpeter@601
   147
    int _arc_num;
kpeter@601
   148
kpeter@601
   149
    // Node and arc maps
kpeter@601
   150
    CapacityVector _cap;
kpeter@601
   151
    CostVector _cost;
kpeter@601
   152
    CostVector _supply;
kpeter@601
   153
    CapacityVector _flow;
kpeter@601
   154
    CostVector _pi;
kpeter@601
   155
kpeter@601
   156
    // Node and arc maps for the spanning tree structure
kpeter@601
   157
    IntVector _depth;
kpeter@601
   158
    IntVector _parent;
kpeter@601
   159
    IntVector _pred;
kpeter@601
   160
    IntVector _thread;
kpeter@601
   161
    BoolVector _forward;
kpeter@601
   162
    IntVector _state;
kpeter@601
   163
kpeter@601
   164
    // The root node
kpeter@601
   165
    int _root;
kpeter@601
   166
kpeter@601
   167
    // The entering arc in the current pivot iteration
kpeter@601
   168
    int _in_arc;
kpeter@601
   169
kpeter@601
   170
    // Temporary data used in the current pivot iteration
kpeter@601
   171
    int join, u_in, v_in, u_out, v_out;
kpeter@601
   172
    int right, first, second, last;
kpeter@601
   173
    int stem, par_stem, new_stem;
kpeter@601
   174
    Capacity delta;
kpeter@601
   175
kpeter@601
   176
  private:
kpeter@601
   177
kpeter@601
   178
    /// \brief Implementation of the "First Eligible" pivot rule for the
kpeter@601
   179
    /// \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   180
    ///
kpeter@601
   181
    /// This class implements the "First Eligible" pivot rule
kpeter@601
   182
    /// for the \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   183
    ///
kpeter@601
   184
    /// For more information see \ref NetworkSimplex::run().
kpeter@601
   185
    class FirstEligiblePivotRule
kpeter@601
   186
    {
kpeter@601
   187
    private:
kpeter@601
   188
kpeter@601
   189
      // References to the NetworkSimplex class
kpeter@601
   190
      const ArcVector &_arc;
kpeter@601
   191
      const IntVector  &_source;
kpeter@601
   192
      const IntVector  &_target;
kpeter@601
   193
      const CostVector &_cost;
kpeter@601
   194
      const IntVector  &_state;
kpeter@601
   195
      const CostVector &_pi;
kpeter@601
   196
      int &_in_arc;
kpeter@601
   197
      int _arc_num;
kpeter@601
   198
kpeter@601
   199
      // Pivot rule data
kpeter@601
   200
      int _next_arc;
kpeter@601
   201
kpeter@601
   202
    public:
kpeter@601
   203
kpeter@601
   204
      /// Constructor
kpeter@601
   205
      FirstEligiblePivotRule(NetworkSimplex &ns) :
kpeter@601
   206
        _arc(ns._arc), _source(ns._source), _target(ns._target),
kpeter@601
   207
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@601
   208
        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
kpeter@601
   209
      {}
kpeter@601
   210
kpeter@601
   211
      /// Find next entering arc
kpeter@601
   212
      bool findEnteringArc() {
kpeter@601
   213
        Cost c;
kpeter@601
   214
        for (int e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   215
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   216
          if (c < 0) {
kpeter@601
   217
            _in_arc = e;
kpeter@601
   218
            _next_arc = e + 1;
kpeter@601
   219
            return true;
kpeter@601
   220
          }
kpeter@601
   221
        }
kpeter@601
   222
        for (int e = 0; e < _next_arc; ++e) {
kpeter@601
   223
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   224
          if (c < 0) {
kpeter@601
   225
            _in_arc = e;
kpeter@601
   226
            _next_arc = e + 1;
kpeter@601
   227
            return true;
kpeter@601
   228
          }
kpeter@601
   229
        }
kpeter@601
   230
        return false;
kpeter@601
   231
      }
kpeter@601
   232
kpeter@601
   233
    }; //class FirstEligiblePivotRule
kpeter@601
   234
kpeter@601
   235
kpeter@601
   236
    /// \brief Implementation of the "Best Eligible" pivot rule for the
kpeter@601
   237
    /// \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   238
    ///
kpeter@601
   239
    /// This class implements the "Best Eligible" pivot rule
kpeter@601
   240
    /// for the \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   241
    ///
kpeter@601
   242
    /// For more information see \ref NetworkSimplex::run().
kpeter@601
   243
    class BestEligiblePivotRule
kpeter@601
   244
    {
kpeter@601
   245
    private:
kpeter@601
   246
kpeter@601
   247
      // References to the NetworkSimplex class
kpeter@601
   248
      const ArcVector &_arc;
kpeter@601
   249
      const IntVector  &_source;
kpeter@601
   250
      const IntVector  &_target;
kpeter@601
   251
      const CostVector &_cost;
kpeter@601
   252
      const IntVector  &_state;
kpeter@601
   253
      const CostVector &_pi;
kpeter@601
   254
      int &_in_arc;
kpeter@601
   255
      int _arc_num;
kpeter@601
   256
kpeter@601
   257
    public:
kpeter@601
   258
kpeter@601
   259
      /// Constructor
kpeter@601
   260
      BestEligiblePivotRule(NetworkSimplex &ns) :
kpeter@601
   261
        _arc(ns._arc), _source(ns._source), _target(ns._target),
kpeter@601
   262
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@601
   263
        _in_arc(ns._in_arc), _arc_num(ns._arc_num)
kpeter@601
   264
      {}
kpeter@601
   265
kpeter@601
   266
      /// Find next entering arc
kpeter@601
   267
      bool findEnteringArc() {
kpeter@601
   268
        Cost c, min = 0;
kpeter@601
   269
        for (int e = 0; e < _arc_num; ++e) {
kpeter@601
   270
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   271
          if (c < min) {
kpeter@601
   272
            min = c;
kpeter@601
   273
            _in_arc = e;
kpeter@601
   274
          }
kpeter@601
   275
        }
kpeter@601
   276
        return min < 0;
kpeter@601
   277
      }
kpeter@601
   278
kpeter@601
   279
    }; //class BestEligiblePivotRule
kpeter@601
   280
kpeter@601
   281
kpeter@601
   282
    /// \brief Implementation of the "Block Search" pivot rule for the
kpeter@601
   283
    /// \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   284
    ///
kpeter@601
   285
    /// This class implements the "Block Search" pivot rule
kpeter@601
   286
    /// for the \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   287
    ///
kpeter@601
   288
    /// For more information see \ref NetworkSimplex::run().
kpeter@601
   289
    class BlockSearchPivotRule
kpeter@601
   290
    {
kpeter@601
   291
    private:
kpeter@601
   292
kpeter@601
   293
      // References to the NetworkSimplex class
kpeter@601
   294
      const ArcVector &_arc;
kpeter@601
   295
      const IntVector  &_source;
kpeter@601
   296
      const IntVector  &_target;
kpeter@601
   297
      const CostVector &_cost;
kpeter@601
   298
      const IntVector  &_state;
kpeter@601
   299
      const CostVector &_pi;
kpeter@601
   300
      int &_in_arc;
kpeter@601
   301
      int _arc_num;
kpeter@601
   302
kpeter@601
   303
      // Pivot rule data
kpeter@601
   304
      int _block_size;
kpeter@601
   305
      int _next_arc;
kpeter@601
   306
kpeter@601
   307
    public:
kpeter@601
   308
kpeter@601
   309
      /// Constructor
kpeter@601
   310
      BlockSearchPivotRule(NetworkSimplex &ns) :
kpeter@601
   311
        _arc(ns._arc), _source(ns._source), _target(ns._target),
kpeter@601
   312
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@601
   313
        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
kpeter@601
   314
      {
kpeter@601
   315
        // The main parameters of the pivot rule
kpeter@601
   316
        const double BLOCK_SIZE_FACTOR = 2.0;
kpeter@601
   317
        const int MIN_BLOCK_SIZE = 10;
kpeter@601
   318
kpeter@601
   319
        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
kpeter@601
   320
                                MIN_BLOCK_SIZE );
kpeter@601
   321
      }
kpeter@601
   322
kpeter@601
   323
      /// Find next entering arc
kpeter@601
   324
      bool findEnteringArc() {
kpeter@601
   325
        Cost c, min = 0;
kpeter@601
   326
        int cnt = _block_size;
kpeter@601
   327
        int e, min_arc = _next_arc;
kpeter@601
   328
        for (e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   329
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   330
          if (c < min) {
kpeter@601
   331
            min = c;
kpeter@601
   332
            min_arc = e;
kpeter@601
   333
          }
kpeter@601
   334
          if (--cnt == 0) {
kpeter@601
   335
            if (min < 0) break;
kpeter@601
   336
            cnt = _block_size;
kpeter@601
   337
          }
kpeter@601
   338
        }
kpeter@601
   339
        if (min == 0 || cnt > 0) {
kpeter@601
   340
          for (e = 0; e < _next_arc; ++e) {
kpeter@601
   341
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   342
            if (c < min) {
kpeter@601
   343
              min = c;
kpeter@601
   344
              min_arc = e;
kpeter@601
   345
            }
kpeter@601
   346
            if (--cnt == 0) {
kpeter@601
   347
              if (min < 0) break;
kpeter@601
   348
              cnt = _block_size;
kpeter@601
   349
            }
kpeter@601
   350
          }
kpeter@601
   351
        }
kpeter@601
   352
        if (min >= 0) return false;
kpeter@601
   353
        _in_arc = min_arc;
kpeter@601
   354
        _next_arc = e;
kpeter@601
   355
        return true;
kpeter@601
   356
      }
kpeter@601
   357
kpeter@601
   358
    }; //class BlockSearchPivotRule
kpeter@601
   359
kpeter@601
   360
kpeter@601
   361
    /// \brief Implementation of the "Candidate List" pivot rule for the
kpeter@601
   362
    /// \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   363
    ///
kpeter@601
   364
    /// This class implements the "Candidate List" pivot rule
kpeter@601
   365
    /// for the \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   366
    ///
kpeter@601
   367
    /// For more information see \ref NetworkSimplex::run().
kpeter@601
   368
    class CandidateListPivotRule
kpeter@601
   369
    {
kpeter@601
   370
    private:
kpeter@601
   371
kpeter@601
   372
      // References to the NetworkSimplex class
kpeter@601
   373
      const ArcVector &_arc;
kpeter@601
   374
      const IntVector  &_source;
kpeter@601
   375
      const IntVector  &_target;
kpeter@601
   376
      const CostVector &_cost;
kpeter@601
   377
      const IntVector  &_state;
kpeter@601
   378
      const CostVector &_pi;
kpeter@601
   379
      int &_in_arc;
kpeter@601
   380
      int _arc_num;
kpeter@601
   381
kpeter@601
   382
      // Pivot rule data
kpeter@601
   383
      IntVector _candidates;
kpeter@601
   384
      int _list_length, _minor_limit;
kpeter@601
   385
      int _curr_length, _minor_count;
kpeter@601
   386
      int _next_arc;
kpeter@601
   387
kpeter@601
   388
    public:
kpeter@601
   389
kpeter@601
   390
      /// Constructor
kpeter@601
   391
      CandidateListPivotRule(NetworkSimplex &ns) :
kpeter@601
   392
        _arc(ns._arc), _source(ns._source), _target(ns._target),
kpeter@601
   393
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@601
   394
        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
kpeter@601
   395
      {
kpeter@601
   396
        // The main parameters of the pivot rule
kpeter@601
   397
        const double LIST_LENGTH_FACTOR = 1.0;
kpeter@601
   398
        const int MIN_LIST_LENGTH = 10;
kpeter@601
   399
        const double MINOR_LIMIT_FACTOR = 0.1;
kpeter@601
   400
        const int MIN_MINOR_LIMIT = 3;
kpeter@601
   401
kpeter@601
   402
        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
kpeter@601
   403
                                 MIN_LIST_LENGTH );
kpeter@601
   404
        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
kpeter@601
   405
                                 MIN_MINOR_LIMIT );
kpeter@601
   406
        _curr_length = _minor_count = 0;
kpeter@601
   407
        _candidates.resize(_list_length);
kpeter@601
   408
      }
kpeter@601
   409
kpeter@601
   410
      /// Find next entering arc
kpeter@601
   411
      bool findEnteringArc() {
kpeter@601
   412
        Cost min, c;
kpeter@601
   413
        int e, min_arc = _next_arc;
kpeter@601
   414
        if (_curr_length > 0 && _minor_count < _minor_limit) {
kpeter@601
   415
          // Minor iteration: select the best eligible arc from the
kpeter@601
   416
          // current candidate list
kpeter@601
   417
          ++_minor_count;
kpeter@601
   418
          min = 0;
kpeter@601
   419
          for (int i = 0; i < _curr_length; ++i) {
kpeter@601
   420
            e = _candidates[i];
kpeter@601
   421
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   422
            if (c < min) {
kpeter@601
   423
              min = c;
kpeter@601
   424
              min_arc = e;
kpeter@601
   425
            }
kpeter@601
   426
            if (c >= 0) {
kpeter@601
   427
              _candidates[i--] = _candidates[--_curr_length];
kpeter@601
   428
            }
kpeter@601
   429
          }
kpeter@601
   430
          if (min < 0) {
kpeter@601
   431
            _in_arc = min_arc;
kpeter@601
   432
            return true;
kpeter@601
   433
          }
kpeter@601
   434
        }
kpeter@601
   435
kpeter@601
   436
        // Major iteration: build a new candidate list
kpeter@601
   437
        min = 0;
kpeter@601
   438
        _curr_length = 0;
kpeter@601
   439
        for (e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   440
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   441
          if (c < 0) {
kpeter@601
   442
            _candidates[_curr_length++] = e;
kpeter@601
   443
            if (c < min) {
kpeter@601
   444
              min = c;
kpeter@601
   445
              min_arc = e;
kpeter@601
   446
            }
kpeter@601
   447
            if (_curr_length == _list_length) break;
kpeter@601
   448
          }
kpeter@601
   449
        }
kpeter@601
   450
        if (_curr_length < _list_length) {
kpeter@601
   451
          for (e = 0; e < _next_arc; ++e) {
kpeter@601
   452
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   453
            if (c < 0) {
kpeter@601
   454
              _candidates[_curr_length++] = e;
kpeter@601
   455
              if (c < min) {
kpeter@601
   456
                min = c;
kpeter@601
   457
                min_arc = e;
kpeter@601
   458
              }
kpeter@601
   459
              if (_curr_length == _list_length) break;
kpeter@601
   460
            }
kpeter@601
   461
          }
kpeter@601
   462
        }
kpeter@601
   463
        if (_curr_length == 0) return false;
kpeter@601
   464
        _minor_count = 1;
kpeter@601
   465
        _in_arc = min_arc;
kpeter@601
   466
        _next_arc = e;
kpeter@601
   467
        return true;
kpeter@601
   468
      }
kpeter@601
   469
kpeter@601
   470
    }; //class CandidateListPivotRule
kpeter@601
   471
kpeter@601
   472
kpeter@601
   473
    /// \brief Implementation of the "Altering Candidate List" pivot rule
kpeter@601
   474
    /// for the \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   475
    ///
kpeter@601
   476
    /// This class implements the "Altering Candidate List" pivot rule
kpeter@601
   477
    /// for the \ref NetworkSimplex "network simplex" algorithm.
kpeter@601
   478
    ///
kpeter@601
   479
    /// For more information see \ref NetworkSimplex::run().
kpeter@601
   480
    class AlteringListPivotRule
kpeter@601
   481
    {
kpeter@601
   482
    private:
kpeter@601
   483
kpeter@601
   484
      // References to the NetworkSimplex class
kpeter@601
   485
      const ArcVector &_arc;
kpeter@601
   486
      const IntVector  &_source;
kpeter@601
   487
      const IntVector  &_target;
kpeter@601
   488
      const CostVector &_cost;
kpeter@601
   489
      const IntVector  &_state;
kpeter@601
   490
      const CostVector &_pi;
kpeter@601
   491
      int &_in_arc;
kpeter@601
   492
      int _arc_num;
kpeter@601
   493
kpeter@601
   494
      // Pivot rule data
kpeter@601
   495
      int _block_size, _head_length, _curr_length;
kpeter@601
   496
      int _next_arc;
kpeter@601
   497
      IntVector _candidates;
kpeter@601
   498
      CostVector _cand_cost;
kpeter@601
   499
kpeter@601
   500
      // Functor class to compare arcs during sort of the candidate list
kpeter@601
   501
      class SortFunc
kpeter@601
   502
      {
kpeter@601
   503
      private:
kpeter@601
   504
        const CostVector &_map;
kpeter@601
   505
      public:
kpeter@601
   506
        SortFunc(const CostVector &map) : _map(map) {}
kpeter@601
   507
        bool operator()(int left, int right) {
kpeter@601
   508
          return _map[left] > _map[right];
kpeter@601
   509
        }
kpeter@601
   510
      };
kpeter@601
   511
kpeter@601
   512
      SortFunc _sort_func;
kpeter@601
   513
kpeter@601
   514
    public:
kpeter@601
   515
kpeter@601
   516
      /// Constructor
kpeter@601
   517
      AlteringListPivotRule(NetworkSimplex &ns) :
kpeter@601
   518
        _arc(ns._arc), _source(ns._source), _target(ns._target),
kpeter@601
   519
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@601
   520
        _in_arc(ns._in_arc), _arc_num(ns._arc_num),
kpeter@601
   521
        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
kpeter@601
   522
      {
kpeter@601
   523
        // The main parameters of the pivot rule
kpeter@601
   524
        const double BLOCK_SIZE_FACTOR = 1.5;
kpeter@601
   525
        const int MIN_BLOCK_SIZE = 10;
kpeter@601
   526
        const double HEAD_LENGTH_FACTOR = 0.1;
kpeter@601
   527
        const int MIN_HEAD_LENGTH = 3;
kpeter@601
   528
kpeter@601
   529
        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
kpeter@601
   530
                                MIN_BLOCK_SIZE );
kpeter@601
   531
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
kpeter@601
   532
                                 MIN_HEAD_LENGTH );
kpeter@601
   533
        _candidates.resize(_head_length + _block_size);
kpeter@601
   534
        _curr_length = 0;
kpeter@601
   535
      }
kpeter@601
   536
kpeter@601
   537
      /// Find next entering arc
kpeter@601
   538
      bool findEnteringArc() {
kpeter@601
   539
        // Check the current candidate list
kpeter@601
   540
        int e;
kpeter@601
   541
        for (int i = 0; i < _curr_length; ++i) {
kpeter@601
   542
          e = _candidates[i];
kpeter@601
   543
          _cand_cost[e] = _state[e] *
kpeter@601
   544
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   545
          if (_cand_cost[e] >= 0) {
kpeter@601
   546
            _candidates[i--] = _candidates[--_curr_length];
kpeter@601
   547
          }
kpeter@601
   548
        }
kpeter@601
   549
kpeter@601
   550
        // Extend the list
kpeter@601
   551
        int cnt = _block_size;
kpeter@601
   552
        int last_edge = 0;
kpeter@601
   553
        int limit = _head_length;
kpeter@601
   554
kpeter@601
   555
        for (int e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   556
          _cand_cost[e] = _state[e] *
kpeter@601
   557
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   558
          if (_cand_cost[e] < 0) {
kpeter@601
   559
            _candidates[_curr_length++] = e;
kpeter@601
   560
            last_edge = e;
kpeter@601
   561
          }
kpeter@601
   562
          if (--cnt == 0) {
kpeter@601
   563
            if (_curr_length > limit) break;
kpeter@601
   564
            limit = 0;
kpeter@601
   565
            cnt = _block_size;
kpeter@601
   566
          }
kpeter@601
   567
        }
kpeter@601
   568
        if (_curr_length <= limit) {
kpeter@601
   569
          for (int e = 0; e < _next_arc; ++e) {
kpeter@601
   570
            _cand_cost[e] = _state[e] *
kpeter@601
   571
              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   572
            if (_cand_cost[e] < 0) {
kpeter@601
   573
              _candidates[_curr_length++] = e;
kpeter@601
   574
              last_edge = e;
kpeter@601
   575
            }
kpeter@601
   576
            if (--cnt == 0) {
kpeter@601
   577
              if (_curr_length > limit) break;
kpeter@601
   578
              limit = 0;
kpeter@601
   579
              cnt = _block_size;
kpeter@601
   580
            }
kpeter@601
   581
          }
kpeter@601
   582
        }
kpeter@601
   583
        if (_curr_length == 0) return false;
kpeter@601
   584
        _next_arc = last_edge + 1;
kpeter@601
   585
kpeter@601
   586
        // Make heap of the candidate list (approximating a partial sort)
kpeter@601
   587
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
kpeter@601
   588
                   _sort_func );
kpeter@601
   589
kpeter@601
   590
        // Pop the first element of the heap
kpeter@601
   591
        _in_arc = _candidates[0];
kpeter@601
   592
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
kpeter@601
   593
                  _sort_func );
kpeter@601
   594
        _curr_length = std::min(_head_length, _curr_length - 1);
kpeter@601
   595
        return true;
kpeter@601
   596
      }
kpeter@601
   597
kpeter@601
   598
    }; //class AlteringListPivotRule
kpeter@601
   599
kpeter@601
   600
  public:
kpeter@601
   601
kpeter@601
   602
    /// \brief General constructor (with lower bounds).
kpeter@601
   603
    ///
kpeter@601
   604
    /// General constructor (with lower bounds).
kpeter@601
   605
    ///
kpeter@601
   606
    /// \param digraph The digraph the algorithm runs on.
kpeter@601
   607
    /// \param lower The lower bounds of the arcs.
kpeter@601
   608
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@601
   609
    /// \param cost The cost (length) values of the arcs.
kpeter@601
   610
    /// \param supply The supply values of the nodes (signed).
kpeter@601
   611
    NetworkSimplex( const Digraph &digraph,
kpeter@601
   612
                    const LowerMap &lower,
kpeter@601
   613
                    const CapacityMap &capacity,
kpeter@601
   614
                    const CostMap &cost,
kpeter@601
   615
                    const SupplyMap &supply ) :
kpeter@601
   616
      _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity),
kpeter@601
   617
      _orig_cost(cost), _orig_supply(&supply),
kpeter@601
   618
      _flow_result(NULL), _potential_result(NULL),
kpeter@601
   619
      _local_flow(false), _local_potential(false),
kpeter@601
   620
      _node_id(digraph)
kpeter@601
   621
    {}
kpeter@601
   622
kpeter@601
   623
    /// \brief General constructor (without lower bounds).
kpeter@601
   624
    ///
kpeter@601
   625
    /// General constructor (without lower bounds).
kpeter@601
   626
    ///
kpeter@601
   627
    /// \param digraph The digraph the algorithm runs on.
kpeter@601
   628
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@601
   629
    /// \param cost The cost (length) values of the arcs.
kpeter@601
   630
    /// \param supply The supply values of the nodes (signed).
kpeter@601
   631
    NetworkSimplex( const Digraph &digraph,
kpeter@601
   632
                    const CapacityMap &capacity,
kpeter@601
   633
                    const CostMap &cost,
kpeter@601
   634
                    const SupplyMap &supply ) :
kpeter@601
   635
      _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity),
kpeter@601
   636
      _orig_cost(cost), _orig_supply(&supply),
kpeter@601
   637
      _flow_result(NULL), _potential_result(NULL),
kpeter@601
   638
      _local_flow(false), _local_potential(false),
kpeter@601
   639
      _node_id(digraph)
kpeter@601
   640
    {}
kpeter@601
   641
kpeter@601
   642
    /// \brief Simple constructor (with lower bounds).
kpeter@601
   643
    ///
kpeter@601
   644
    /// Simple constructor (with lower bounds).
kpeter@601
   645
    ///
kpeter@601
   646
    /// \param digraph The digraph the algorithm runs on.
kpeter@601
   647
    /// \param lower The lower bounds of the arcs.
kpeter@601
   648
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@601
   649
    /// \param cost The cost (length) values of the arcs.
kpeter@601
   650
    /// \param s The source node.
kpeter@601
   651
    /// \param t The target node.
kpeter@601
   652
    /// \param flow_value The required amount of flow from node \c s
kpeter@601
   653
    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
kpeter@601
   654
    NetworkSimplex( const Digraph &digraph,
kpeter@601
   655
                    const LowerMap &lower,
kpeter@601
   656
                    const CapacityMap &capacity,
kpeter@601
   657
                    const CostMap &cost,
kpeter@601
   658
                    Node s, Node t,
kpeter@601
   659
                    Capacity flow_value ) :
kpeter@601
   660
      _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity),
kpeter@601
   661
      _orig_cost(cost), _orig_supply(NULL),
kpeter@601
   662
      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
kpeter@601
   663
      _flow_result(NULL), _potential_result(NULL),
kpeter@601
   664
      _local_flow(false), _local_potential(false),
kpeter@601
   665
      _node_id(digraph)
kpeter@601
   666
    {}
kpeter@601
   667
kpeter@601
   668
    /// \brief Simple constructor (without lower bounds).
kpeter@601
   669
    ///
kpeter@601
   670
    /// Simple constructor (without lower bounds).
kpeter@601
   671
    ///
kpeter@601
   672
    /// \param digraph The digraph the algorithm runs on.
kpeter@601
   673
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@601
   674
    /// \param cost The cost (length) values of the arcs.
kpeter@601
   675
    /// \param s The source node.
kpeter@601
   676
    /// \param t The target node.
kpeter@601
   677
    /// \param flow_value The required amount of flow from node \c s
kpeter@601
   678
    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
kpeter@601
   679
    NetworkSimplex( const Digraph &digraph,
kpeter@601
   680
                    const CapacityMap &capacity,
kpeter@601
   681
                    const CostMap &cost,
kpeter@601
   682
                    Node s, Node t,
kpeter@601
   683
                    Capacity flow_value ) :
kpeter@601
   684
      _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity),
kpeter@601
   685
      _orig_cost(cost), _orig_supply(NULL),
kpeter@601
   686
      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
kpeter@601
   687
      _flow_result(NULL), _potential_result(NULL),
kpeter@601
   688
      _local_flow(false), _local_potential(false),
kpeter@601
   689
      _node_id(digraph)
kpeter@601
   690
    {}
kpeter@601
   691
kpeter@601
   692
    /// Destructor.
kpeter@601
   693
    ~NetworkSimplex() {
kpeter@601
   694
      if (_local_flow) delete _flow_result;
kpeter@601
   695
      if (_local_potential) delete _potential_result;
kpeter@601
   696
    }
kpeter@601
   697
kpeter@601
   698
    /// \brief Set the flow map.
kpeter@601
   699
    ///
kpeter@601
   700
    /// This function sets the flow map.
kpeter@601
   701
    ///
kpeter@601
   702
    /// \return <tt>(*this)</tt>
kpeter@601
   703
    NetworkSimplex& flowMap(FlowMap &map) {
kpeter@601
   704
      if (_local_flow) {
kpeter@601
   705
        delete _flow_result;
kpeter@601
   706
        _local_flow = false;
kpeter@601
   707
      }
kpeter@601
   708
      _flow_result = &map;
kpeter@601
   709
      return *this;
kpeter@601
   710
    }
kpeter@601
   711
kpeter@601
   712
    /// \brief Set the potential map.
kpeter@601
   713
    ///
kpeter@601
   714
    /// This function sets the potential map.
kpeter@601
   715
    ///
kpeter@601
   716
    /// \return <tt>(*this)</tt>
kpeter@601
   717
    NetworkSimplex& potentialMap(PotentialMap &map) {
kpeter@601
   718
      if (_local_potential) {
kpeter@601
   719
        delete _potential_result;
kpeter@601
   720
        _local_potential = false;
kpeter@601
   721
      }
kpeter@601
   722
      _potential_result = &map;
kpeter@601
   723
      return *this;
kpeter@601
   724
    }
kpeter@601
   725
kpeter@601
   726
    /// \name Execution control
kpeter@601
   727
    /// The algorithm can be executed using the
kpeter@601
   728
    /// \ref lemon::NetworkSimplex::run() "run()" function.
kpeter@601
   729
    /// @{
kpeter@601
   730
kpeter@601
   731
    /// \brief Run the algorithm.
kpeter@601
   732
    ///
kpeter@601
   733
    /// This function runs the algorithm.
kpeter@601
   734
    ///
kpeter@601
   735
    /// \param pivot_rule The pivot rule that is used during the
kpeter@601
   736
    /// algorithm.
kpeter@601
   737
    ///
kpeter@601
   738
    /// The available pivot rules:
kpeter@601
   739
    ///
kpeter@601
   740
    /// - FIRST_ELIGIBLE_PIVOT The next eligible arc is selected in
kpeter@601
   741
    /// a wraparound fashion in every iteration
kpeter@601
   742
    /// (\ref FirstEligiblePivotRule).
kpeter@601
   743
    ///
kpeter@601
   744
    /// - BEST_ELIGIBLE_PIVOT The best eligible arc is selected in
kpeter@601
   745
    /// every iteration (\ref BestEligiblePivotRule).
kpeter@601
   746
    ///
kpeter@601
   747
    /// - BLOCK_SEARCH_PIVOT A specified number of arcs are examined in
kpeter@601
   748
    /// every iteration in a wraparound fashion and the best eligible
kpeter@601
   749
    /// arc is selected from this block (\ref BlockSearchPivotRule).
kpeter@601
   750
    ///
kpeter@601
   751
    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
kpeter@601
   752
    /// built from eligible arcs in a wraparound fashion and in the
kpeter@601
   753
    /// following minor iterations the best eligible arc is selected
kpeter@601
   754
    /// from this list (\ref CandidateListPivotRule).
kpeter@601
   755
    ///
kpeter@601
   756
    /// - ALTERING_LIST_PIVOT It is a modified version of the
kpeter@601
   757
    /// "Candidate List" pivot rule. It keeps only the several best
kpeter@601
   758
    /// eligible arcs from the former candidate list and extends this
kpeter@601
   759
    /// list in every iteration (\ref AlteringListPivotRule).
kpeter@601
   760
    ///
kpeter@601
   761
    /// According to our comprehensive benchmark tests the "Block Search"
kpeter@601
   762
    /// pivot rule proved to be the fastest and the most robust on
kpeter@601
   763
    /// various test inputs. Thus it is the default option.
kpeter@601
   764
    ///
kpeter@601
   765
    /// \return \c true if a feasible flow can be found.
kpeter@601
   766
    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
kpeter@601
   767
      return init() && start(pivot_rule);
kpeter@601
   768
    }
kpeter@601
   769
kpeter@601
   770
    /// @}
kpeter@601
   771
kpeter@601
   772
    /// \name Query Functions
kpeter@601
   773
    /// The results of the algorithm can be obtained using these
kpeter@601
   774
    /// functions.\n
kpeter@601
   775
    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
kpeter@601
   776
    /// using them.
kpeter@601
   777
    /// @{
kpeter@601
   778
kpeter@601
   779
    /// \brief Return a const reference to the flow map.
kpeter@601
   780
    ///
kpeter@601
   781
    /// This function returns a const reference to an arc map storing
kpeter@601
   782
    /// the found flow.
kpeter@601
   783
    ///
kpeter@601
   784
    /// \pre \ref run() must be called before using this function.
kpeter@601
   785
    const FlowMap& flowMap() const {
kpeter@601
   786
      return *_flow_result;
kpeter@601
   787
    }
kpeter@601
   788
kpeter@601
   789
    /// \brief Return a const reference to the potential map
kpeter@601
   790
    /// (the dual solution).
kpeter@601
   791
    ///
kpeter@601
   792
    /// This function returns a const reference to a node map storing
kpeter@601
   793
    /// the found potentials (the dual solution).
kpeter@601
   794
    ///
kpeter@601
   795
    /// \pre \ref run() must be called before using this function.
kpeter@601
   796
    const PotentialMap& potentialMap() const {
kpeter@601
   797
      return *_potential_result;
kpeter@601
   798
    }
kpeter@601
   799
kpeter@601
   800
    /// \brief Return the flow on the given arc.
kpeter@601
   801
    ///
kpeter@601
   802
    /// This function returns the flow on the given arc.
kpeter@601
   803
    ///
kpeter@601
   804
    /// \pre \ref run() must be called before using this function.
kpeter@601
   805
    Capacity flow(const Arc& arc) const {
kpeter@601
   806
      return (*_flow_result)[arc];
kpeter@601
   807
    }
kpeter@601
   808
kpeter@601
   809
    /// \brief Return the potential of the given node.
kpeter@601
   810
    ///
kpeter@601
   811
    /// This function returns the potential of the given node.
kpeter@601
   812
    ///
kpeter@601
   813
    /// \pre \ref run() must be called before using this function.
kpeter@601
   814
    Cost potential(const Node& node) const {
kpeter@601
   815
      return (*_potential_result)[node];
kpeter@601
   816
    }
kpeter@601
   817
kpeter@601
   818
    /// \brief Return the total cost of the found flow.
kpeter@601
   819
    ///
kpeter@601
   820
    /// This function returns the total cost of the found flow.
kpeter@601
   821
    /// The complexity of the function is \f$ O(e) \f$.
kpeter@601
   822
    ///
kpeter@601
   823
    /// \pre \ref run() must be called before using this function.
kpeter@601
   824
    Cost totalCost() const {
kpeter@601
   825
      Cost c = 0;
kpeter@601
   826
      for (ArcIt e(_orig_graph); e != INVALID; ++e)
kpeter@601
   827
        c += (*_flow_result)[e] * _orig_cost[e];
kpeter@601
   828
      return c;
kpeter@601
   829
    }
kpeter@601
   830
kpeter@601
   831
    /// @}
kpeter@601
   832
kpeter@601
   833
  private:
kpeter@601
   834
kpeter@601
   835
    // Initialize internal data structures
kpeter@601
   836
    bool init() {
kpeter@601
   837
      // Initialize result maps
kpeter@601
   838
      if (!_flow_result) {
kpeter@601
   839
        _flow_result = new FlowMap(_orig_graph);
kpeter@601
   840
        _local_flow = true;
kpeter@601
   841
      }
kpeter@601
   842
      if (!_potential_result) {
kpeter@601
   843
        _potential_result = new PotentialMap(_orig_graph);
kpeter@601
   844
        _local_potential = true;
kpeter@601
   845
      }
kpeter@601
   846
kpeter@601
   847
      // Initialize vectors
kpeter@601
   848
      _node_num = countNodes(_orig_graph);
kpeter@601
   849
      _arc_num = countArcs(_orig_graph);
kpeter@601
   850
      int all_node_num = _node_num + 1;
kpeter@601
   851
      int all_edge_num = _arc_num + _node_num;
kpeter@601
   852
kpeter@601
   853
      _arc.resize(_arc_num);
kpeter@601
   854
      _node.reserve(_node_num);
kpeter@601
   855
      _source.resize(all_edge_num);
kpeter@601
   856
      _target.resize(all_edge_num);
kpeter@601
   857
kpeter@601
   858
      _cap.resize(all_edge_num);
kpeter@601
   859
      _cost.resize(all_edge_num);
kpeter@601
   860
      _supply.resize(all_node_num);
kpeter@601
   861
      _flow.resize(all_edge_num, 0);
kpeter@601
   862
      _pi.resize(all_node_num, 0);
kpeter@601
   863
kpeter@601
   864
      _depth.resize(all_node_num);
kpeter@601
   865
      _parent.resize(all_node_num);
kpeter@601
   866
      _pred.resize(all_node_num);
kpeter@601
   867
      _thread.resize(all_node_num);
kpeter@601
   868
      _forward.resize(all_node_num);
kpeter@601
   869
      _state.resize(all_edge_num, STATE_LOWER);
kpeter@601
   870
kpeter@601
   871
      // Initialize node related data
kpeter@601
   872
      bool valid_supply = true;
kpeter@601
   873
      if (_orig_supply) {
kpeter@601
   874
        Supply sum = 0;
kpeter@601
   875
        int i = 0;
kpeter@601
   876
        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
kpeter@601
   877
          _node.push_back(n);
kpeter@601
   878
          _node_id[n] = i;
kpeter@601
   879
          _supply[i] = (*_orig_supply)[n];
kpeter@601
   880
          sum += _supply[i];
kpeter@601
   881
        }
kpeter@601
   882
        valid_supply = (sum == 0);
kpeter@601
   883
      } else {
kpeter@601
   884
        int i = 0;
kpeter@601
   885
        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
kpeter@601
   886
          _node.push_back(n);
kpeter@601
   887
          _node_id[n] = i;
kpeter@601
   888
          _supply[i] = 0;
kpeter@601
   889
        }
kpeter@601
   890
        _supply[_node_id[_orig_source]] =  _orig_flow_value;
kpeter@601
   891
        _supply[_node_id[_orig_target]] = -_orig_flow_value;
kpeter@601
   892
      }
kpeter@601
   893
      if (!valid_supply) return false;
kpeter@601
   894
kpeter@601
   895
      // Set data for the artificial root node
kpeter@601
   896
      _root = _node_num;
kpeter@601
   897
      _depth[_root] = 0;
kpeter@601
   898
      _parent[_root] = -1;
kpeter@601
   899
      _pred[_root] = -1;
kpeter@601
   900
      _thread[_root] = 0;
kpeter@601
   901
      _supply[_root] = 0;
kpeter@601
   902
      _pi[_root] = 0;
kpeter@601
   903
kpeter@601
   904
      // Store the arcs in a mixed order
kpeter@601
   905
      int k = std::max(int(sqrt(_arc_num)), 10);
kpeter@601
   906
      int i = 0;
kpeter@601
   907
      for (ArcIt e(_orig_graph); e != INVALID; ++e) {
kpeter@601
   908
        _arc[i] = e;
kpeter@601
   909
        if ((i += k) >= _arc_num) i = (i % k) + 1;
kpeter@601
   910
      }
kpeter@601
   911
kpeter@601
   912
      // Initialize arc maps
kpeter@601
   913
      for (int i = 0; i != _arc_num; ++i) {
kpeter@601
   914
        Arc e = _arc[i];
kpeter@601
   915
        _source[i] = _node_id[_orig_graph.source(e)];
kpeter@601
   916
        _target[i] = _node_id[_orig_graph.target(e)];
kpeter@601
   917
        _cost[i] = _orig_cost[e];
kpeter@601
   918
        _cap[i] = _orig_cap[e];
kpeter@601
   919
      }
kpeter@601
   920
kpeter@601
   921
      // Remove non-zero lower bounds
kpeter@601
   922
      if (_orig_lower) {
kpeter@601
   923
        for (int i = 0; i != _arc_num; ++i) {
kpeter@601
   924
          Capacity c = (*_orig_lower)[_arc[i]];
kpeter@601
   925
          if (c != 0) {
kpeter@601
   926
            _cap[i] -= c;
kpeter@601
   927
            _supply[_source[i]] -= c;
kpeter@601
   928
            _supply[_target[i]] += c;
kpeter@601
   929
          }
kpeter@601
   930
        }
kpeter@601
   931
      }
kpeter@601
   932
kpeter@601
   933
      // Add artificial arcs and initialize the spanning tree data structure
kpeter@601
   934
      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
kpeter@601
   935
      Capacity max_cap = std::numeric_limits<Capacity>::max();
kpeter@601
   936
      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
kpeter@601
   937
        _thread[u] = u + 1;
kpeter@601
   938
        _depth[u] = 1;
kpeter@601
   939
        _parent[u] = _root;
kpeter@601
   940
        _pred[u] = e;
kpeter@601
   941
        if (_supply[u] >= 0) {
kpeter@601
   942
          _flow[e] = _supply[u];
kpeter@601
   943
          _forward[u] = true;
kpeter@601
   944
          _pi[u] = -max_cost;
kpeter@601
   945
        } else {
kpeter@601
   946
          _flow[e] = -_supply[u];
kpeter@601
   947
          _forward[u] = false;
kpeter@601
   948
          _pi[u] = max_cost;
kpeter@601
   949
        }
kpeter@601
   950
        _cost[e] = max_cost;
kpeter@601
   951
        _cap[e] = max_cap;
kpeter@601
   952
        _state[e] = STATE_TREE;
kpeter@601
   953
      }
kpeter@601
   954
kpeter@601
   955
      return true;
kpeter@601
   956
    }
kpeter@601
   957
kpeter@601
   958
    // Find the join node
kpeter@601
   959
    void findJoinNode() {
kpeter@601
   960
      int u = _source[_in_arc];
kpeter@601
   961
      int v = _target[_in_arc];
kpeter@601
   962
      while (_depth[u] > _depth[v]) u = _parent[u];
kpeter@601
   963
      while (_depth[v] > _depth[u]) v = _parent[v];
kpeter@601
   964
      while (u != v) {
kpeter@601
   965
        u = _parent[u];
kpeter@601
   966
        v = _parent[v];
kpeter@601
   967
      }
kpeter@601
   968
      join = u;
kpeter@601
   969
    }
kpeter@601
   970
kpeter@601
   971
    // Find the leaving arc of the cycle and returns true if the
kpeter@601
   972
    // leaving arc is not the same as the entering arc
kpeter@601
   973
    bool findLeavingArc() {
kpeter@601
   974
      // Initialize first and second nodes according to the direction
kpeter@601
   975
      // of the cycle
kpeter@601
   976
      if (_state[_in_arc] == STATE_LOWER) {
kpeter@601
   977
        first  = _source[_in_arc];
kpeter@601
   978
        second = _target[_in_arc];
kpeter@601
   979
      } else {
kpeter@601
   980
        first  = _target[_in_arc];
kpeter@601
   981
        second = _source[_in_arc];
kpeter@601
   982
      }
kpeter@601
   983
      delta = _cap[_in_arc];
kpeter@601
   984
      int result = 0;
kpeter@601
   985
      Capacity d;
kpeter@601
   986
      int e;
kpeter@601
   987
kpeter@601
   988
      // Search the cycle along the path form the first node to the root
kpeter@601
   989
      for (int u = first; u != join; u = _parent[u]) {
kpeter@601
   990
        e = _pred[u];
kpeter@601
   991
        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
kpeter@601
   992
        if (d < delta) {
kpeter@601
   993
          delta = d;
kpeter@601
   994
          u_out = u;
kpeter@601
   995
          result = 1;
kpeter@601
   996
        }
kpeter@601
   997
      }
kpeter@601
   998
      // Search the cycle along the path form the second node to the root
kpeter@601
   999
      for (int u = second; u != join; u = _parent[u]) {
kpeter@601
  1000
        e = _pred[u];
kpeter@601
  1001
        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
kpeter@601
  1002
        if (d <= delta) {
kpeter@601
  1003
          delta = d;
kpeter@601
  1004
          u_out = u;
kpeter@601
  1005
          result = 2;
kpeter@601
  1006
        }
kpeter@601
  1007
      }
kpeter@601
  1008
kpeter@601
  1009
      if (result == 1) {
kpeter@601
  1010
        u_in = first;
kpeter@601
  1011
        v_in = second;
kpeter@601
  1012
      } else {
kpeter@601
  1013
        u_in = second;
kpeter@601
  1014
        v_in = first;
kpeter@601
  1015
      }
kpeter@601
  1016
      return result != 0;
kpeter@601
  1017
    }
kpeter@601
  1018
kpeter@601
  1019
    // Change _flow and _state vectors
kpeter@601
  1020
    void changeFlow(bool change) {
kpeter@601
  1021
      // Augment along the cycle
kpeter@601
  1022
      if (delta > 0) {
kpeter@601
  1023
        Capacity val = _state[_in_arc] * delta;
kpeter@601
  1024
        _flow[_in_arc] += val;
kpeter@601
  1025
        for (int u = _source[_in_arc]; u != join; u = _parent[u]) {
kpeter@601
  1026
          _flow[_pred[u]] += _forward[u] ? -val : val;
kpeter@601
  1027
        }
kpeter@601
  1028
        for (int u = _target[_in_arc]; u != join; u = _parent[u]) {
kpeter@601
  1029
          _flow[_pred[u]] += _forward[u] ? val : -val;
kpeter@601
  1030
        }
kpeter@601
  1031
      }
kpeter@601
  1032
      // Update the state of the entering and leaving arcs
kpeter@601
  1033
      if (change) {
kpeter@601
  1034
        _state[_in_arc] = STATE_TREE;
kpeter@601
  1035
        _state[_pred[u_out]] =
kpeter@601
  1036
          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
kpeter@601
  1037
      } else {
kpeter@601
  1038
        _state[_in_arc] = -_state[_in_arc];
kpeter@601
  1039
      }
kpeter@601
  1040
    }
kpeter@601
  1041
kpeter@601
  1042
    // Update _thread and _parent vectors
kpeter@601
  1043
    void updateThreadParent() {
kpeter@601
  1044
      int u;
kpeter@601
  1045
      v_out = _parent[u_out];
kpeter@601
  1046
kpeter@601
  1047
      // Handle the case when join and v_out coincide
kpeter@601
  1048
      bool par_first = false;
kpeter@601
  1049
      if (join == v_out) {
kpeter@601
  1050
        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
kpeter@601
  1051
        if (u == v_in) {
kpeter@601
  1052
          par_first = true;
kpeter@601
  1053
          while (_thread[u] != u_out) u = _thread[u];
kpeter@601
  1054
          first = u;
kpeter@601
  1055
        }
kpeter@601
  1056
      }
kpeter@601
  1057
kpeter@601
  1058
      // Find the last successor of u_in (u) and the node after it (right)
kpeter@601
  1059
      // according to the thread index
kpeter@601
  1060
      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
kpeter@601
  1061
      right = _thread[u];
kpeter@601
  1062
      if (_thread[v_in] == u_out) {
kpeter@601
  1063
        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
kpeter@601
  1064
        if (last == u_out) last = _thread[last];
kpeter@601
  1065
      }
kpeter@601
  1066
      else last = _thread[v_in];
kpeter@601
  1067
kpeter@601
  1068
      // Update stem nodes
kpeter@601
  1069
      _thread[v_in] = stem = u_in;
kpeter@601
  1070
      par_stem = v_in;
kpeter@601
  1071
      while (stem != u_out) {
kpeter@601
  1072
        _thread[u] = new_stem = _parent[stem];
kpeter@601
  1073
kpeter@601
  1074
        // Find the node just before the stem node (u) according to
kpeter@601
  1075
        // the original thread index
kpeter@601
  1076
        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
kpeter@601
  1077
        _thread[u] = right;
kpeter@601
  1078
kpeter@601
  1079
        // Change the parent node of stem and shift stem and par_stem nodes
kpeter@601
  1080
        _parent[stem] = par_stem;
kpeter@601
  1081
        par_stem = stem;
kpeter@601
  1082
        stem = new_stem;
kpeter@601
  1083
kpeter@601
  1084
        // Find the last successor of stem (u) and the node after it (right)
kpeter@601
  1085
        // according to the thread index
kpeter@601
  1086
        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
kpeter@601
  1087
        right = _thread[u];
kpeter@601
  1088
      }
kpeter@601
  1089
      _parent[u_out] = par_stem;
kpeter@601
  1090
      _thread[u] = last;
kpeter@601
  1091
kpeter@601
  1092
      if (join == v_out && par_first) {
kpeter@601
  1093
        if (first != v_in) _thread[first] = right;
kpeter@601
  1094
      } else {
kpeter@601
  1095
        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
kpeter@601
  1096
        _thread[u] = right;
kpeter@601
  1097
      }
kpeter@601
  1098
    }
kpeter@601
  1099
kpeter@601
  1100
    // Update _pred and _forward vectors
kpeter@601
  1101
    void updatePredArc() {
kpeter@601
  1102
      int u = u_out, v;
kpeter@601
  1103
      while (u != u_in) {
kpeter@601
  1104
        v = _parent[u];
kpeter@601
  1105
        _pred[u] = _pred[v];
kpeter@601
  1106
        _forward[u] = !_forward[v];
kpeter@601
  1107
        u = v;
kpeter@601
  1108
      }
kpeter@601
  1109
      _pred[u_in] = _in_arc;
kpeter@601
  1110
      _forward[u_in] = (u_in == _source[_in_arc]);
kpeter@601
  1111
    }
kpeter@601
  1112
kpeter@601
  1113
    // Update _depth and _potential vectors
kpeter@601
  1114
    void updateDepthPotential() {
kpeter@601
  1115
      _depth[u_in] = _depth[v_in] + 1;
kpeter@601
  1116
      Cost sigma = _forward[u_in] ?
kpeter@601
  1117
        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
kpeter@601
  1118
        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
kpeter@601
  1119
      _pi[u_in] += sigma;
kpeter@601
  1120
      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
kpeter@601
  1121
        _depth[u] = _depth[_parent[u]] + 1;
kpeter@601
  1122
        if (_depth[u] <= _depth[u_in]) break;
kpeter@601
  1123
        _pi[u] += sigma;
kpeter@601
  1124
      }
kpeter@601
  1125
    }
kpeter@601
  1126
kpeter@601
  1127
    // Execute the algorithm
kpeter@601
  1128
    bool start(PivotRuleEnum pivot_rule) {
kpeter@601
  1129
      // Select the pivot rule implementation
kpeter@601
  1130
      switch (pivot_rule) {
kpeter@601
  1131
        case FIRST_ELIGIBLE_PIVOT:
kpeter@601
  1132
          return start<FirstEligiblePivotRule>();
kpeter@601
  1133
        case BEST_ELIGIBLE_PIVOT:
kpeter@601
  1134
          return start<BestEligiblePivotRule>();
kpeter@601
  1135
        case BLOCK_SEARCH_PIVOT:
kpeter@601
  1136
          return start<BlockSearchPivotRule>();
kpeter@601
  1137
        case CANDIDATE_LIST_PIVOT:
kpeter@601
  1138
          return start<CandidateListPivotRule>();
kpeter@601
  1139
        case ALTERING_LIST_PIVOT:
kpeter@601
  1140
          return start<AlteringListPivotRule>();
kpeter@601
  1141
      }
kpeter@601
  1142
      return false;
kpeter@601
  1143
    }
kpeter@601
  1144
kpeter@601
  1145
    template<class PivotRuleImplementation>
kpeter@601
  1146
    bool start() {
kpeter@601
  1147
      PivotRuleImplementation pivot(*this);
kpeter@601
  1148
kpeter@601
  1149
      // Execute the network simplex algorithm
kpeter@601
  1150
      while (pivot.findEnteringArc()) {
kpeter@601
  1151
        findJoinNode();
kpeter@601
  1152
        bool change = findLeavingArc();
kpeter@601
  1153
        changeFlow(change);
kpeter@601
  1154
        if (change) {
kpeter@601
  1155
          updateThreadParent();
kpeter@601
  1156
          updatePredArc();
kpeter@601
  1157
          updateDepthPotential();
kpeter@601
  1158
        }
kpeter@601
  1159
      }
kpeter@601
  1160
kpeter@601
  1161
      // Check if the flow amount equals zero on all the artificial arcs
kpeter@601
  1162
      for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
kpeter@601
  1163
        if (_flow[e] > 0) return false;
kpeter@601
  1164
      }
kpeter@601
  1165
kpeter@601
  1166
      // Copy flow values to _flow_result
kpeter@601
  1167
      if (_orig_lower) {
kpeter@601
  1168
        for (int i = 0; i != _arc_num; ++i) {
kpeter@601
  1169
          Arc e = _arc[i];
kpeter@601
  1170
          (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
kpeter@601
  1171
        }
kpeter@601
  1172
      } else {
kpeter@601
  1173
        for (int i = 0; i != _arc_num; ++i) {
kpeter@601
  1174
          (*_flow_result)[_arc[i]] = _flow[i];
kpeter@601
  1175
        }
kpeter@601
  1176
      }
kpeter@601
  1177
      // Copy potential values to _potential_result
kpeter@601
  1178
      for (int i = 0; i != _node_num; ++i) {
kpeter@601
  1179
        (*_potential_result)[_node[i]] = _pi[i];
kpeter@601
  1180
      }
kpeter@601
  1181
kpeter@601
  1182
      return true;
kpeter@601
  1183
    }
kpeter@601
  1184
kpeter@601
  1185
  }; //class NetworkSimplex
kpeter@601
  1186
kpeter@601
  1187
  ///@}
kpeter@601
  1188
kpeter@601
  1189
} //namespace lemon
kpeter@601
  1190
kpeter@601
  1191
#endif //LEMON_NETWORK_SIMPLEX_H