lemon/network_simplex.h
author Peter Kovacs <kpeter@inf.elte.hu>
Tue, 24 Mar 2009 00:18:25 +0100
changeset 651 8c3112a66878
parent 650 425cc8328c0e
child 652 5232721b3f14
permissions -rw-r--r--
Use XTI implementation instead of ATI in NetworkSimplex (#234)

XTI (eXtended Threaded Index) is an imporved version of the widely
known ATI (Augmented Threaded Index) method for storing and updating
the spanning tree structure in Network Simplex algorithms.

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