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