lemon/network_simplex.h
author Akos Ladanyi <ladanyi@tmit.bme.hu>
Thu, 23 Apr 2009 07:30:40 +0100
changeset 621 b536eaacb39b
parent 608 6ac5d9ae1d3d
child 612 0c8e5c688440
child 613 b1811c363299
permissions -rw-r--r--
FindCOIN for CMake (#256)
kpeter@601
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
kpeter@601
     2
 *
kpeter@601
     3
 * This file is a part of LEMON, a generic C++ optimization library.
kpeter@601
     4
 *
kpeter@601
     5
 * Copyright (C) 2003-2009
kpeter@601
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
kpeter@601
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
kpeter@601
     8
 *
kpeter@601
     9
 * Permission to use, modify and distribute this software is granted
kpeter@601
    10
 * provided that this copyright notice appears in all copies. For
kpeter@601
    11
 * precise terms see the accompanying LICENSE file.
kpeter@601
    12
 *
kpeter@601
    13
 * This software is provided "AS IS" with no warranty of any kind,
kpeter@601
    14
 * express or implied, and with no claim as to its suitability for any
kpeter@601
    15
 * purpose.
kpeter@601
    16
 *
kpeter@601
    17
 */
kpeter@601
    18
kpeter@601
    19
#ifndef LEMON_NETWORK_SIMPLEX_H
kpeter@601
    20
#define LEMON_NETWORK_SIMPLEX_H
kpeter@601
    21
kpeter@601
    22
/// \ingroup min_cost_flow
kpeter@601
    23
///
kpeter@601
    24
/// \file
kpeter@605
    25
/// \brief Network Simplex algorithm for finding a minimum cost flow.
kpeter@601
    26
kpeter@601
    27
#include <vector>
kpeter@601
    28
#include <limits>
kpeter@601
    29
#include <algorithm>
kpeter@601
    30
kpeter@603
    31
#include <lemon/core.h>
kpeter@601
    32
#include <lemon/math.h>
kpeter@609
    33
#include <lemon/maps.h>
kpeter@609
    34
#include <lemon/circulation.h>
kpeter@609
    35
#include <lemon/adaptors.h>
kpeter@601
    36
kpeter@601
    37
namespace lemon {
kpeter@601
    38
kpeter@601
    39
  /// \addtogroup min_cost_flow
kpeter@601
    40
  /// @{
kpeter@601
    41
kpeter@605
    42
  /// \brief Implementation of the primal Network Simplex algorithm
kpeter@601
    43
  /// for finding a \ref min_cost_flow "minimum cost flow".
kpeter@601
    44
  ///
kpeter@605
    45
  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
kpeter@601
    46
  /// for finding a \ref min_cost_flow "minimum cost flow".
kpeter@606
    47
  /// This algorithm is a specialized version of the linear programming
kpeter@606
    48
  /// simplex method directly for the minimum cost flow problem.
kpeter@606
    49
  /// It is one of the most efficient solution methods.
kpeter@606
    50
  ///
kpeter@606
    51
  /// In general this class is the fastest implementation available
kpeter@606
    52
  /// in LEMON for the minimum cost flow problem.
kpeter@609
    53
  /// Moreover it supports both direction of the supply/demand inequality
kpeter@609
    54
  /// constraints. For more information see \ref ProblemType.
kpeter@601
    55
  ///
kpeter@605
    56
  /// \tparam GR The digraph type the algorithm runs on.
kpeter@607
    57
  /// \tparam F The value type used for flow amounts, capacity bounds
kpeter@607
    58
  /// and supply values in the algorithm. By default it is \c int.
kpeter@607
    59
  /// \tparam C The value type used for costs and potentials in the
kpeter@607
    60
  /// algorithm. By default it is the same as \c F.
kpeter@601
    61
  ///
kpeter@608
    62
  /// \warning Both value types must be signed and all input data must
kpeter@608
    63
  /// be integer.
kpeter@601
    64
  ///
kpeter@605
    65
  /// \note %NetworkSimplex provides five different pivot rule
kpeter@609
    66
  /// implementations, from which the most efficient one is used
kpeter@609
    67
  /// by default. For more information see \ref PivotRule.
kpeter@607
    68
  template <typename GR, typename F = int, typename C = F>
kpeter@601
    69
  class NetworkSimplex
kpeter@601
    70
  {
kpeter@605
    71
  public:
kpeter@601
    72
kpeter@607
    73
    /// The flow type of the algorithm
kpeter@607
    74
    typedef F Flow;
kpeter@607
    75
    /// The cost type of the algorithm
kpeter@607
    76
    typedef C Cost;
kpeter@609
    77
#ifdef DOXYGEN
kpeter@609
    78
    /// The type of the flow map
kpeter@609
    79
    typedef GR::ArcMap<Flow> FlowMap;
kpeter@609
    80
    /// The type of the potential map
kpeter@609
    81
    typedef GR::NodeMap<Cost> PotentialMap;
kpeter@609
    82
#else
kpeter@605
    83
    /// The type of the flow map
kpeter@607
    84
    typedef typename GR::template ArcMap<Flow> FlowMap;
kpeter@605
    85
    /// The type of the potential map
kpeter@607
    86
    typedef typename GR::template NodeMap<Cost> PotentialMap;
kpeter@609
    87
#endif
kpeter@605
    88
kpeter@605
    89
  public:
kpeter@605
    90
kpeter@605
    91
    /// \brief Enum type for selecting the pivot rule.
kpeter@605
    92
    ///
kpeter@605
    93
    /// Enum type for selecting the pivot rule for the \ref run()
kpeter@605
    94
    /// function.
kpeter@605
    95
    ///
kpeter@605
    96
    /// \ref NetworkSimplex provides five different pivot rule
kpeter@605
    97
    /// implementations that significantly affect the running time
kpeter@605
    98
    /// of the algorithm.
kpeter@605
    99
    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
kpeter@605
   100
    /// proved to be the most efficient and the most robust on various
kpeter@605
   101
    /// test inputs according to our benchmark tests.
kpeter@605
   102
    /// However another pivot rule can be selected using the \ref run()
kpeter@605
   103
    /// function with the proper parameter.
kpeter@605
   104
    enum PivotRule {
kpeter@605
   105
kpeter@605
   106
      /// The First Eligible pivot rule.
kpeter@605
   107
      /// The next eligible arc is selected in a wraparound fashion
kpeter@605
   108
      /// in every iteration.
kpeter@605
   109
      FIRST_ELIGIBLE,
kpeter@605
   110
kpeter@605
   111
      /// The Best Eligible pivot rule.
kpeter@605
   112
      /// The best eligible arc is selected in every iteration.
kpeter@605
   113
      BEST_ELIGIBLE,
kpeter@605
   114
kpeter@605
   115
      /// The Block Search pivot rule.
kpeter@605
   116
      /// A specified number of arcs are examined in every iteration
kpeter@605
   117
      /// in a wraparound fashion and the best eligible arc is selected
kpeter@605
   118
      /// from this block.
kpeter@605
   119
      BLOCK_SEARCH,
kpeter@605
   120
kpeter@605
   121
      /// The Candidate List pivot rule.
kpeter@605
   122
      /// In a major iteration a candidate list is built from eligible arcs
kpeter@605
   123
      /// in a wraparound fashion and in the following minor iterations
kpeter@605
   124
      /// the best eligible arc is selected from this list.
kpeter@605
   125
      CANDIDATE_LIST,
kpeter@605
   126
kpeter@605
   127
      /// The Altering Candidate List pivot rule.
kpeter@605
   128
      /// It is a modified version of the Candidate List method.
kpeter@605
   129
      /// It keeps only the several best eligible arcs from the former
kpeter@605
   130
      /// candidate list and extends this list in every iteration.
kpeter@605
   131
      ALTERING_LIST
kpeter@605
   132
    };
kpeter@609
   133
    
kpeter@609
   134
    /// \brief Enum type for selecting the problem type.
kpeter@609
   135
    ///
kpeter@609
   136
    /// Enum type for selecting the problem type, i.e. the direction of
kpeter@609
   137
    /// the inequalities in the supply/demand constraints of the
kpeter@609
   138
    /// \ref min_cost_flow "minimum cost flow problem".
kpeter@609
   139
    ///
kpeter@609
   140
    /// The default problem type is \c GEQ, since this form is supported
kpeter@609
   141
    /// by other minimum cost flow algorithms and the \ref Circulation
kpeter@609
   142
    /// algorithm as well.
kpeter@609
   143
    /// The \c LEQ problem type can be selected using the \ref problemType()
kpeter@609
   144
    /// function.
kpeter@609
   145
    ///
kpeter@609
   146
    /// Note that the equality form is a special case of both problem type.
kpeter@609
   147
    enum ProblemType {
kpeter@609
   148
kpeter@609
   149
      /// This option means that there are "<em>greater or equal</em>"
kpeter@609
   150
      /// constraints in the defintion, i.e. the exact formulation of the
kpeter@609
   151
      /// problem is the following.
kpeter@609
   152
      /**
kpeter@609
   153
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
kpeter@609
   154
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
kpeter@609
   155
              sup(u) \quad \forall u\in V \f]
kpeter@609
   156
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
kpeter@609
   157
      */
kpeter@609
   158
      /// It means that the total demand must be greater or equal to the 
kpeter@609
   159
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
kpeter@609
   160
      /// negative) and all the supplies have to be carried out from 
kpeter@609
   161
      /// the supply nodes, but there could be demands that are not 
kpeter@609
   162
      /// satisfied.
kpeter@609
   163
      GEQ,
kpeter@609
   164
      /// It is just an alias for the \c GEQ option.
kpeter@609
   165
      CARRY_SUPPLIES = GEQ,
kpeter@609
   166
kpeter@609
   167
      /// This option means that there are "<em>less or equal</em>"
kpeter@609
   168
      /// constraints in the defintion, i.e. the exact formulation of the
kpeter@609
   169
      /// problem is the following.
kpeter@609
   170
      /**
kpeter@609
   171
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
kpeter@609
   172
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
kpeter@609
   173
              sup(u) \quad \forall u\in V \f]
kpeter@609
   174
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
kpeter@609
   175
      */
kpeter@609
   176
      /// It means that the total demand must be less or equal to the 
kpeter@609
   177
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
kpeter@609
   178
      /// positive) and all the demands have to be satisfied, but there
kpeter@609
   179
      /// could be supplies that are not carried out from the supply
kpeter@609
   180
      /// nodes.
kpeter@609
   181
      LEQ,
kpeter@609
   182
      /// It is just an alias for the \c LEQ option.
kpeter@609
   183
      SATISFY_DEMANDS = LEQ
kpeter@609
   184
    };
kpeter@605
   185
kpeter@605
   186
  private:
kpeter@605
   187
kpeter@605
   188
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
kpeter@605
   189
kpeter@607
   190
    typedef typename GR::template ArcMap<Flow> FlowArcMap;
kpeter@607
   191
    typedef typename GR::template ArcMap<Cost> CostArcMap;
kpeter@607
   192
    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
kpeter@601
   193
kpeter@601
   194
    typedef std::vector<Arc> ArcVector;
kpeter@601
   195
    typedef std::vector<Node> NodeVector;
kpeter@601
   196
    typedef std::vector<int> IntVector;
kpeter@601
   197
    typedef std::vector<bool> BoolVector;
kpeter@607
   198
    typedef std::vector<Flow> FlowVector;
kpeter@607
   199
    typedef std::vector<Cost> CostVector;
kpeter@601
   200
kpeter@601
   201
    // State constants for arcs
kpeter@601
   202
    enum ArcStateEnum {
kpeter@601
   203
      STATE_UPPER = -1,
kpeter@601
   204
      STATE_TREE  =  0,
kpeter@601
   205
      STATE_LOWER =  1
kpeter@601
   206
    };
kpeter@601
   207
kpeter@601
   208
  private:
kpeter@601
   209
kpeter@605
   210
    // Data related to the underlying digraph
kpeter@605
   211
    const GR &_graph;
kpeter@605
   212
    int _node_num;
kpeter@605
   213
    int _arc_num;
kpeter@605
   214
kpeter@605
   215
    // Parameters of the problem
kpeter@607
   216
    FlowArcMap *_plower;
kpeter@607
   217
    FlowArcMap *_pupper;
kpeter@607
   218
    CostArcMap *_pcost;
kpeter@607
   219
    FlowNodeMap *_psupply;
kpeter@605
   220
    bool _pstsup;
kpeter@605
   221
    Node _psource, _ptarget;
kpeter@607
   222
    Flow _pstflow;
kpeter@609
   223
    ProblemType _ptype;
kpeter@601
   224
kpeter@601
   225
    // Result maps
kpeter@603
   226
    FlowMap *_flow_map;
kpeter@603
   227
    PotentialMap *_potential_map;
kpeter@601
   228
    bool _local_flow;
kpeter@601
   229
    bool _local_potential;
kpeter@601
   230
kpeter@605
   231
    // Data structures for storing the digraph
kpeter@603
   232
    IntNodeMap _node_id;
kpeter@603
   233
    ArcVector _arc_ref;
kpeter@603
   234
    IntVector _source;
kpeter@603
   235
    IntVector _target;
kpeter@603
   236
kpeter@605
   237
    // Node and arc data
kpeter@607
   238
    FlowVector _cap;
kpeter@607
   239
    CostVector _cost;
kpeter@607
   240
    FlowVector _supply;
kpeter@607
   241
    FlowVector _flow;
kpeter@607
   242
    CostVector _pi;
kpeter@601
   243
kpeter@603
   244
    // Data for storing the spanning tree structure
kpeter@601
   245
    IntVector _parent;
kpeter@601
   246
    IntVector _pred;
kpeter@601
   247
    IntVector _thread;
kpeter@604
   248
    IntVector _rev_thread;
kpeter@604
   249
    IntVector _succ_num;
kpeter@604
   250
    IntVector _last_succ;
kpeter@604
   251
    IntVector _dirty_revs;
kpeter@601
   252
    BoolVector _forward;
kpeter@601
   253
    IntVector _state;
kpeter@601
   254
    int _root;
kpeter@601
   255
kpeter@601
   256
    // Temporary data used in the current pivot iteration
kpeter@603
   257
    int in_arc, join, u_in, v_in, u_out, v_out;
kpeter@603
   258
    int first, second, right, last;
kpeter@601
   259
    int stem, par_stem, new_stem;
kpeter@607
   260
    Flow delta;
kpeter@601
   261
kpeter@601
   262
  private:
kpeter@601
   263
kpeter@605
   264
    // Implementation of the First Eligible pivot rule
kpeter@601
   265
    class FirstEligiblePivotRule
kpeter@601
   266
    {
kpeter@601
   267
    private:
kpeter@601
   268
kpeter@601
   269
      // References to the NetworkSimplex class
kpeter@601
   270
      const IntVector  &_source;
kpeter@601
   271
      const IntVector  &_target;
kpeter@607
   272
      const CostVector &_cost;
kpeter@601
   273
      const IntVector  &_state;
kpeter@607
   274
      const CostVector &_pi;
kpeter@601
   275
      int &_in_arc;
kpeter@601
   276
      int _arc_num;
kpeter@601
   277
kpeter@601
   278
      // Pivot rule data
kpeter@601
   279
      int _next_arc;
kpeter@601
   280
kpeter@601
   281
    public:
kpeter@601
   282
kpeter@605
   283
      // Constructor
kpeter@601
   284
      FirstEligiblePivotRule(NetworkSimplex &ns) :
kpeter@603
   285
        _source(ns._source), _target(ns._target),
kpeter@601
   286
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@603
   287
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
kpeter@601
   288
      {}
kpeter@601
   289
kpeter@605
   290
      // Find next entering arc
kpeter@601
   291
      bool findEnteringArc() {
kpeter@607
   292
        Cost c;
kpeter@601
   293
        for (int e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   294
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   295
          if (c < 0) {
kpeter@601
   296
            _in_arc = e;
kpeter@601
   297
            _next_arc = e + 1;
kpeter@601
   298
            return true;
kpeter@601
   299
          }
kpeter@601
   300
        }
kpeter@601
   301
        for (int e = 0; e < _next_arc; ++e) {
kpeter@601
   302
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   303
          if (c < 0) {
kpeter@601
   304
            _in_arc = e;
kpeter@601
   305
            _next_arc = e + 1;
kpeter@601
   306
            return true;
kpeter@601
   307
          }
kpeter@601
   308
        }
kpeter@601
   309
        return false;
kpeter@601
   310
      }
kpeter@601
   311
kpeter@601
   312
    }; //class FirstEligiblePivotRule
kpeter@601
   313
kpeter@601
   314
kpeter@605
   315
    // Implementation of the Best Eligible pivot rule
kpeter@601
   316
    class BestEligiblePivotRule
kpeter@601
   317
    {
kpeter@601
   318
    private:
kpeter@601
   319
kpeter@601
   320
      // References to the NetworkSimplex class
kpeter@601
   321
      const IntVector  &_source;
kpeter@601
   322
      const IntVector  &_target;
kpeter@607
   323
      const CostVector &_cost;
kpeter@601
   324
      const IntVector  &_state;
kpeter@607
   325
      const CostVector &_pi;
kpeter@601
   326
      int &_in_arc;
kpeter@601
   327
      int _arc_num;
kpeter@601
   328
kpeter@601
   329
    public:
kpeter@601
   330
kpeter@605
   331
      // Constructor
kpeter@601
   332
      BestEligiblePivotRule(NetworkSimplex &ns) :
kpeter@603
   333
        _source(ns._source), _target(ns._target),
kpeter@601
   334
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@603
   335
        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
kpeter@601
   336
      {}
kpeter@601
   337
kpeter@605
   338
      // Find next entering arc
kpeter@601
   339
      bool findEnteringArc() {
kpeter@607
   340
        Cost c, min = 0;
kpeter@601
   341
        for (int e = 0; e < _arc_num; ++e) {
kpeter@601
   342
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   343
          if (c < min) {
kpeter@601
   344
            min = c;
kpeter@601
   345
            _in_arc = e;
kpeter@601
   346
          }
kpeter@601
   347
        }
kpeter@601
   348
        return min < 0;
kpeter@601
   349
      }
kpeter@601
   350
kpeter@601
   351
    }; //class BestEligiblePivotRule
kpeter@601
   352
kpeter@601
   353
kpeter@605
   354
    // Implementation of the Block Search pivot rule
kpeter@601
   355
    class BlockSearchPivotRule
kpeter@601
   356
    {
kpeter@601
   357
    private:
kpeter@601
   358
kpeter@601
   359
      // References to the NetworkSimplex class
kpeter@601
   360
      const IntVector  &_source;
kpeter@601
   361
      const IntVector  &_target;
kpeter@607
   362
      const CostVector &_cost;
kpeter@601
   363
      const IntVector  &_state;
kpeter@607
   364
      const CostVector &_pi;
kpeter@601
   365
      int &_in_arc;
kpeter@601
   366
      int _arc_num;
kpeter@601
   367
kpeter@601
   368
      // Pivot rule data
kpeter@601
   369
      int _block_size;
kpeter@601
   370
      int _next_arc;
kpeter@601
   371
kpeter@601
   372
    public:
kpeter@601
   373
kpeter@605
   374
      // Constructor
kpeter@601
   375
      BlockSearchPivotRule(NetworkSimplex &ns) :
kpeter@603
   376
        _source(ns._source), _target(ns._target),
kpeter@601
   377
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@603
   378
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
kpeter@601
   379
      {
kpeter@601
   380
        // The main parameters of the pivot rule
kpeter@601
   381
        const double BLOCK_SIZE_FACTOR = 2.0;
kpeter@601
   382
        const int MIN_BLOCK_SIZE = 10;
kpeter@601
   383
kpeter@601
   384
        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
kpeter@601
   385
                                MIN_BLOCK_SIZE );
kpeter@601
   386
      }
kpeter@601
   387
kpeter@605
   388
      // Find next entering arc
kpeter@601
   389
      bool findEnteringArc() {
kpeter@607
   390
        Cost c, min = 0;
kpeter@601
   391
        int cnt = _block_size;
kpeter@601
   392
        int e, min_arc = _next_arc;
kpeter@601
   393
        for (e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   394
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   395
          if (c < min) {
kpeter@601
   396
            min = c;
kpeter@601
   397
            min_arc = e;
kpeter@601
   398
          }
kpeter@601
   399
          if (--cnt == 0) {
kpeter@601
   400
            if (min < 0) break;
kpeter@601
   401
            cnt = _block_size;
kpeter@601
   402
          }
kpeter@601
   403
        }
kpeter@601
   404
        if (min == 0 || cnt > 0) {
kpeter@601
   405
          for (e = 0; e < _next_arc; ++e) {
kpeter@601
   406
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   407
            if (c < min) {
kpeter@601
   408
              min = c;
kpeter@601
   409
              min_arc = e;
kpeter@601
   410
            }
kpeter@601
   411
            if (--cnt == 0) {
kpeter@601
   412
              if (min < 0) break;
kpeter@601
   413
              cnt = _block_size;
kpeter@601
   414
            }
kpeter@601
   415
          }
kpeter@601
   416
        }
kpeter@601
   417
        if (min >= 0) return false;
kpeter@601
   418
        _in_arc = min_arc;
kpeter@601
   419
        _next_arc = e;
kpeter@601
   420
        return true;
kpeter@601
   421
      }
kpeter@601
   422
kpeter@601
   423
    }; //class BlockSearchPivotRule
kpeter@601
   424
kpeter@601
   425
kpeter@605
   426
    // Implementation of the Candidate List pivot rule
kpeter@601
   427
    class CandidateListPivotRule
kpeter@601
   428
    {
kpeter@601
   429
    private:
kpeter@601
   430
kpeter@601
   431
      // References to the NetworkSimplex class
kpeter@601
   432
      const IntVector  &_source;
kpeter@601
   433
      const IntVector  &_target;
kpeter@607
   434
      const CostVector &_cost;
kpeter@601
   435
      const IntVector  &_state;
kpeter@607
   436
      const CostVector &_pi;
kpeter@601
   437
      int &_in_arc;
kpeter@601
   438
      int _arc_num;
kpeter@601
   439
kpeter@601
   440
      // Pivot rule data
kpeter@601
   441
      IntVector _candidates;
kpeter@601
   442
      int _list_length, _minor_limit;
kpeter@601
   443
      int _curr_length, _minor_count;
kpeter@601
   444
      int _next_arc;
kpeter@601
   445
kpeter@601
   446
    public:
kpeter@601
   447
kpeter@601
   448
      /// Constructor
kpeter@601
   449
      CandidateListPivotRule(NetworkSimplex &ns) :
kpeter@603
   450
        _source(ns._source), _target(ns._target),
kpeter@601
   451
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@603
   452
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
kpeter@601
   453
      {
kpeter@601
   454
        // The main parameters of the pivot rule
kpeter@601
   455
        const double LIST_LENGTH_FACTOR = 1.0;
kpeter@601
   456
        const int MIN_LIST_LENGTH = 10;
kpeter@601
   457
        const double MINOR_LIMIT_FACTOR = 0.1;
kpeter@601
   458
        const int MIN_MINOR_LIMIT = 3;
kpeter@601
   459
kpeter@601
   460
        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
kpeter@601
   461
                                 MIN_LIST_LENGTH );
kpeter@601
   462
        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
kpeter@601
   463
                                 MIN_MINOR_LIMIT );
kpeter@601
   464
        _curr_length = _minor_count = 0;
kpeter@601
   465
        _candidates.resize(_list_length);
kpeter@601
   466
      }
kpeter@601
   467
kpeter@601
   468
      /// Find next entering arc
kpeter@601
   469
      bool findEnteringArc() {
kpeter@607
   470
        Cost min, c;
kpeter@601
   471
        int e, min_arc = _next_arc;
kpeter@601
   472
        if (_curr_length > 0 && _minor_count < _minor_limit) {
kpeter@601
   473
          // Minor iteration: select the best eligible arc from the
kpeter@601
   474
          // current candidate list
kpeter@601
   475
          ++_minor_count;
kpeter@601
   476
          min = 0;
kpeter@601
   477
          for (int i = 0; i < _curr_length; ++i) {
kpeter@601
   478
            e = _candidates[i];
kpeter@601
   479
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   480
            if (c < min) {
kpeter@601
   481
              min = c;
kpeter@601
   482
              min_arc = e;
kpeter@601
   483
            }
kpeter@601
   484
            if (c >= 0) {
kpeter@601
   485
              _candidates[i--] = _candidates[--_curr_length];
kpeter@601
   486
            }
kpeter@601
   487
          }
kpeter@601
   488
          if (min < 0) {
kpeter@601
   489
            _in_arc = min_arc;
kpeter@601
   490
            return true;
kpeter@601
   491
          }
kpeter@601
   492
        }
kpeter@601
   493
kpeter@601
   494
        // Major iteration: build a new candidate list
kpeter@601
   495
        min = 0;
kpeter@601
   496
        _curr_length = 0;
kpeter@601
   497
        for (e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   498
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   499
          if (c < 0) {
kpeter@601
   500
            _candidates[_curr_length++] = e;
kpeter@601
   501
            if (c < min) {
kpeter@601
   502
              min = c;
kpeter@601
   503
              min_arc = e;
kpeter@601
   504
            }
kpeter@601
   505
            if (_curr_length == _list_length) break;
kpeter@601
   506
          }
kpeter@601
   507
        }
kpeter@601
   508
        if (_curr_length < _list_length) {
kpeter@601
   509
          for (e = 0; e < _next_arc; ++e) {
kpeter@601
   510
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   511
            if (c < 0) {
kpeter@601
   512
              _candidates[_curr_length++] = e;
kpeter@601
   513
              if (c < min) {
kpeter@601
   514
                min = c;
kpeter@601
   515
                min_arc = e;
kpeter@601
   516
              }
kpeter@601
   517
              if (_curr_length == _list_length) break;
kpeter@601
   518
            }
kpeter@601
   519
          }
kpeter@601
   520
        }
kpeter@601
   521
        if (_curr_length == 0) return false;
kpeter@601
   522
        _minor_count = 1;
kpeter@601
   523
        _in_arc = min_arc;
kpeter@601
   524
        _next_arc = e;
kpeter@601
   525
        return true;
kpeter@601
   526
      }
kpeter@601
   527
kpeter@601
   528
    }; //class CandidateListPivotRule
kpeter@601
   529
kpeter@601
   530
kpeter@605
   531
    // Implementation of the Altering Candidate List pivot rule
kpeter@601
   532
    class AlteringListPivotRule
kpeter@601
   533
    {
kpeter@601
   534
    private:
kpeter@601
   535
kpeter@601
   536
      // References to the NetworkSimplex class
kpeter@601
   537
      const IntVector  &_source;
kpeter@601
   538
      const IntVector  &_target;
kpeter@607
   539
      const CostVector &_cost;
kpeter@601
   540
      const IntVector  &_state;
kpeter@607
   541
      const CostVector &_pi;
kpeter@601
   542
      int &_in_arc;
kpeter@601
   543
      int _arc_num;
kpeter@601
   544
kpeter@601
   545
      // Pivot rule data
kpeter@601
   546
      int _block_size, _head_length, _curr_length;
kpeter@601
   547
      int _next_arc;
kpeter@601
   548
      IntVector _candidates;
kpeter@607
   549
      CostVector _cand_cost;
kpeter@601
   550
kpeter@601
   551
      // Functor class to compare arcs during sort of the candidate list
kpeter@601
   552
      class SortFunc
kpeter@601
   553
      {
kpeter@601
   554
      private:
kpeter@607
   555
        const CostVector &_map;
kpeter@601
   556
      public:
kpeter@607
   557
        SortFunc(const CostVector &map) : _map(map) {}
kpeter@601
   558
        bool operator()(int left, int right) {
kpeter@601
   559
          return _map[left] > _map[right];
kpeter@601
   560
        }
kpeter@601
   561
      };
kpeter@601
   562
kpeter@601
   563
      SortFunc _sort_func;
kpeter@601
   564
kpeter@601
   565
    public:
kpeter@601
   566
kpeter@605
   567
      // Constructor
kpeter@601
   568
      AlteringListPivotRule(NetworkSimplex &ns) :
kpeter@603
   569
        _source(ns._source), _target(ns._target),
kpeter@601
   570
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
kpeter@603
   571
        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
kpeter@601
   572
        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
kpeter@601
   573
      {
kpeter@601
   574
        // The main parameters of the pivot rule
kpeter@601
   575
        const double BLOCK_SIZE_FACTOR = 1.5;
kpeter@601
   576
        const int MIN_BLOCK_SIZE = 10;
kpeter@601
   577
        const double HEAD_LENGTH_FACTOR = 0.1;
kpeter@601
   578
        const int MIN_HEAD_LENGTH = 3;
kpeter@601
   579
kpeter@601
   580
        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
kpeter@601
   581
                                MIN_BLOCK_SIZE );
kpeter@601
   582
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
kpeter@601
   583
                                 MIN_HEAD_LENGTH );
kpeter@601
   584
        _candidates.resize(_head_length + _block_size);
kpeter@601
   585
        _curr_length = 0;
kpeter@601
   586
      }
kpeter@601
   587
kpeter@605
   588
      // Find next entering arc
kpeter@601
   589
      bool findEnteringArc() {
kpeter@601
   590
        // Check the current candidate list
kpeter@601
   591
        int e;
kpeter@601
   592
        for (int i = 0; i < _curr_length; ++i) {
kpeter@601
   593
          e = _candidates[i];
kpeter@601
   594
          _cand_cost[e] = _state[e] *
kpeter@601
   595
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   596
          if (_cand_cost[e] >= 0) {
kpeter@601
   597
            _candidates[i--] = _candidates[--_curr_length];
kpeter@601
   598
          }
kpeter@601
   599
        }
kpeter@601
   600
kpeter@601
   601
        // Extend the list
kpeter@601
   602
        int cnt = _block_size;
kpeter@603
   603
        int last_arc = 0;
kpeter@601
   604
        int limit = _head_length;
kpeter@601
   605
kpeter@601
   606
        for (int e = _next_arc; e < _arc_num; ++e) {
kpeter@601
   607
          _cand_cost[e] = _state[e] *
kpeter@601
   608
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   609
          if (_cand_cost[e] < 0) {
kpeter@601
   610
            _candidates[_curr_length++] = e;
kpeter@603
   611
            last_arc = e;
kpeter@601
   612
          }
kpeter@601
   613
          if (--cnt == 0) {
kpeter@601
   614
            if (_curr_length > limit) break;
kpeter@601
   615
            limit = 0;
kpeter@601
   616
            cnt = _block_size;
kpeter@601
   617
          }
kpeter@601
   618
        }
kpeter@601
   619
        if (_curr_length <= limit) {
kpeter@601
   620
          for (int e = 0; e < _next_arc; ++e) {
kpeter@601
   621
            _cand_cost[e] = _state[e] *
kpeter@601
   622
              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
kpeter@601
   623
            if (_cand_cost[e] < 0) {
kpeter@601
   624
              _candidates[_curr_length++] = e;
kpeter@603
   625
              last_arc = e;
kpeter@601
   626
            }
kpeter@601
   627
            if (--cnt == 0) {
kpeter@601
   628
              if (_curr_length > limit) break;
kpeter@601
   629
              limit = 0;
kpeter@601
   630
              cnt = _block_size;
kpeter@601
   631
            }
kpeter@601
   632
          }
kpeter@601
   633
        }
kpeter@601
   634
        if (_curr_length == 0) return false;
kpeter@603
   635
        _next_arc = last_arc + 1;
kpeter@601
   636
kpeter@601
   637
        // Make heap of the candidate list (approximating a partial sort)
kpeter@601
   638
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
kpeter@601
   639
                   _sort_func );
kpeter@601
   640
kpeter@601
   641
        // Pop the first element of the heap
kpeter@601
   642
        _in_arc = _candidates[0];
kpeter@601
   643
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
kpeter@601
   644
                  _sort_func );
kpeter@601
   645
        _curr_length = std::min(_head_length, _curr_length - 1);
kpeter@601
   646
        return true;
kpeter@601
   647
      }
kpeter@601
   648
kpeter@601
   649
    }; //class AlteringListPivotRule
kpeter@601
   650
kpeter@601
   651
  public:
kpeter@601
   652
kpeter@605
   653
    /// \brief Constructor.
kpeter@601
   654
    ///
kpeter@609
   655
    /// The constructor of the class.
kpeter@601
   656
    ///
kpeter@603
   657
    /// \param graph The digraph the algorithm runs on.
kpeter@605
   658
    NetworkSimplex(const GR& graph) :
kpeter@605
   659
      _graph(graph),
kpeter@605
   660
      _plower(NULL), _pupper(NULL), _pcost(NULL),
kpeter@609
   661
      _psupply(NULL), _pstsup(false), _ptype(GEQ),
kpeter@603
   662
      _flow_map(NULL), _potential_map(NULL),
kpeter@601
   663
      _local_flow(false), _local_potential(false),
kpeter@603
   664
      _node_id(graph)
kpeter@605
   665
    {
kpeter@607
   666
      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
kpeter@607
   667
                   std::numeric_limits<Flow>::is_signed,
kpeter@607
   668
        "The flow type of NetworkSimplex must be signed integer");
kpeter@607
   669
      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
kpeter@607
   670
                   std::numeric_limits<Cost>::is_signed,
kpeter@607
   671
        "The cost type of NetworkSimplex must be signed integer");
kpeter@605
   672
    }
kpeter@601
   673
kpeter@601
   674
    /// Destructor.
kpeter@601
   675
    ~NetworkSimplex() {
kpeter@603
   676
      if (_local_flow) delete _flow_map;
kpeter@603
   677
      if (_local_potential) delete _potential_map;
kpeter@601
   678
    }
kpeter@601
   679
kpeter@609
   680
    /// \name Parameters
kpeter@609
   681
    /// The parameters of the algorithm can be specified using these
kpeter@609
   682
    /// functions.
kpeter@609
   683
kpeter@609
   684
    /// @{
kpeter@609
   685
kpeter@605
   686
    /// \brief Set the lower bounds on the arcs.
kpeter@605
   687
    ///
kpeter@605
   688
    /// This function sets the lower bounds on the arcs.
kpeter@605
   689
    /// If neither this function nor \ref boundMaps() is used before
kpeter@605
   690
    /// calling \ref run(), the lower bounds will be set to zero
kpeter@605
   691
    /// on all arcs.
kpeter@605
   692
    ///
kpeter@605
   693
    /// \param map An arc map storing the lower bounds.
kpeter@607
   694
    /// Its \c Value type must be convertible to the \c Flow type
kpeter@605
   695
    /// of the algorithm.
kpeter@605
   696
    ///
kpeter@605
   697
    /// \return <tt>(*this)</tt>
kpeter@605
   698
    template <typename LOWER>
kpeter@605
   699
    NetworkSimplex& lowerMap(const LOWER& map) {
kpeter@605
   700
      delete _plower;
kpeter@607
   701
      _plower = new FlowArcMap(_graph);
kpeter@605
   702
      for (ArcIt a(_graph); a != INVALID; ++a) {
kpeter@605
   703
        (*_plower)[a] = map[a];
kpeter@605
   704
      }
kpeter@605
   705
      return *this;
kpeter@605
   706
    }
kpeter@605
   707
kpeter@605
   708
    /// \brief Set the upper bounds (capacities) on the arcs.
kpeter@605
   709
    ///
kpeter@605
   710
    /// This function sets the upper bounds (capacities) on the arcs.
kpeter@605
   711
    /// If none of the functions \ref upperMap(), \ref capacityMap()
kpeter@605
   712
    /// and \ref boundMaps() is used before calling \ref run(),
kpeter@605
   713
    /// the upper bounds (capacities) will be set to
kpeter@607
   714
    /// \c std::numeric_limits<Flow>::max() on all arcs.
kpeter@605
   715
    ///
kpeter@605
   716
    /// \param map An arc map storing the upper bounds.
kpeter@607
   717
    /// Its \c Value type must be convertible to the \c Flow type
kpeter@605
   718
    /// of the algorithm.
kpeter@605
   719
    ///
kpeter@605
   720
    /// \return <tt>(*this)</tt>
kpeter@605
   721
    template<typename UPPER>
kpeter@605
   722
    NetworkSimplex& upperMap(const UPPER& map) {
kpeter@605
   723
      delete _pupper;
kpeter@607
   724
      _pupper = new FlowArcMap(_graph);
kpeter@605
   725
      for (ArcIt a(_graph); a != INVALID; ++a) {
kpeter@605
   726
        (*_pupper)[a] = map[a];
kpeter@605
   727
      }
kpeter@605
   728
      return *this;
kpeter@605
   729
    }
kpeter@605
   730
kpeter@605
   731
    /// \brief Set the upper bounds (capacities) on the arcs.
kpeter@605
   732
    ///
kpeter@605
   733
    /// This function sets the upper bounds (capacities) on the arcs.
kpeter@605
   734
    /// It is just an alias for \ref upperMap().
kpeter@605
   735
    ///
kpeter@605
   736
    /// \return <tt>(*this)</tt>
kpeter@605
   737
    template<typename CAP>
kpeter@605
   738
    NetworkSimplex& capacityMap(const CAP& map) {
kpeter@605
   739
      return upperMap(map);
kpeter@605
   740
    }
kpeter@605
   741
kpeter@605
   742
    /// \brief Set the lower and upper bounds on the arcs.
kpeter@605
   743
    ///
kpeter@605
   744
    /// This function sets the lower and upper bounds on the arcs.
kpeter@605
   745
    /// If neither this function nor \ref lowerMap() is used before
kpeter@605
   746
    /// calling \ref run(), the lower bounds will be set to zero
kpeter@605
   747
    /// on all arcs.
kpeter@605
   748
    /// If none of the functions \ref upperMap(), \ref capacityMap()
kpeter@605
   749
    /// and \ref boundMaps() is used before calling \ref run(),
kpeter@605
   750
    /// the upper bounds (capacities) will be set to
kpeter@607
   751
    /// \c std::numeric_limits<Flow>::max() on all arcs.
kpeter@605
   752
    ///
kpeter@605
   753
    /// \param lower An arc map storing the lower bounds.
kpeter@605
   754
    /// \param upper An arc map storing the upper bounds.
kpeter@605
   755
    ///
kpeter@605
   756
    /// The \c Value type of the maps must be convertible to the
kpeter@607
   757
    /// \c Flow type of the algorithm.
kpeter@605
   758
    ///
kpeter@605
   759
    /// \note This function is just a shortcut of calling \ref lowerMap()
kpeter@605
   760
    /// and \ref upperMap() separately.
kpeter@605
   761
    ///
kpeter@605
   762
    /// \return <tt>(*this)</tt>
kpeter@605
   763
    template <typename LOWER, typename UPPER>
kpeter@605
   764
    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
kpeter@605
   765
      return lowerMap(lower).upperMap(upper);
kpeter@605
   766
    }
kpeter@605
   767
kpeter@605
   768
    /// \brief Set the costs of the arcs.
kpeter@605
   769
    ///
kpeter@605
   770
    /// This function sets the costs of the arcs.
kpeter@605
   771
    /// If it is not used before calling \ref run(), the costs
kpeter@605
   772
    /// will be set to \c 1 on all arcs.
kpeter@605
   773
    ///
kpeter@605
   774
    /// \param map An arc map storing the costs.
kpeter@607
   775
    /// Its \c Value type must be convertible to the \c Cost type
kpeter@605
   776
    /// of the algorithm.
kpeter@605
   777
    ///
kpeter@605
   778
    /// \return <tt>(*this)</tt>
kpeter@605
   779
    template<typename COST>
kpeter@605
   780
    NetworkSimplex& costMap(const COST& map) {
kpeter@605
   781
      delete _pcost;
kpeter@607
   782
      _pcost = new CostArcMap(_graph);
kpeter@605
   783
      for (ArcIt a(_graph); a != INVALID; ++a) {
kpeter@605
   784
        (*_pcost)[a] = map[a];
kpeter@605
   785
      }
kpeter@605
   786
      return *this;
kpeter@605
   787
    }
kpeter@605
   788
kpeter@605
   789
    /// \brief Set the supply values of the nodes.
kpeter@605
   790
    ///
kpeter@605
   791
    /// This function sets the supply values of the nodes.
kpeter@605
   792
    /// If neither this function nor \ref stSupply() is used before
kpeter@605
   793
    /// calling \ref run(), the supply of each node will be set to zero.
kpeter@605
   794
    /// (It makes sense only if non-zero lower bounds are given.)
kpeter@605
   795
    ///
kpeter@605
   796
    /// \param map A node map storing the supply values.
kpeter@607
   797
    /// Its \c Value type must be convertible to the \c Flow type
kpeter@605
   798
    /// of the algorithm.
kpeter@605
   799
    ///
kpeter@605
   800
    /// \return <tt>(*this)</tt>
kpeter@605
   801
    template<typename SUP>
kpeter@605
   802
    NetworkSimplex& supplyMap(const SUP& map) {
kpeter@605
   803
      delete _psupply;
kpeter@605
   804
      _pstsup = false;
kpeter@607
   805
      _psupply = new FlowNodeMap(_graph);
kpeter@605
   806
      for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@605
   807
        (*_psupply)[n] = map[n];
kpeter@605
   808
      }
kpeter@605
   809
      return *this;
kpeter@605
   810
    }
kpeter@605
   811
kpeter@605
   812
    /// \brief Set single source and target nodes and a supply value.
kpeter@605
   813
    ///
kpeter@605
   814
    /// This function sets a single source node and a single target node
kpeter@605
   815
    /// and the required flow value.
kpeter@605
   816
    /// If neither this function nor \ref supplyMap() is used before
kpeter@605
   817
    /// calling \ref run(), the supply of each node will be set to zero.
kpeter@605
   818
    /// (It makes sense only if non-zero lower bounds are given.)
kpeter@605
   819
    ///
kpeter@605
   820
    /// \param s The source node.
kpeter@605
   821
    /// \param t The target node.
kpeter@605
   822
    /// \param k The required amount of flow from node \c s to node \c t
kpeter@605
   823
    /// (i.e. the supply of \c s and the demand of \c t).
kpeter@605
   824
    ///
kpeter@605
   825
    /// \return <tt>(*this)</tt>
kpeter@607
   826
    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
kpeter@605
   827
      delete _psupply;
kpeter@605
   828
      _psupply = NULL;
kpeter@605
   829
      _pstsup = true;
kpeter@605
   830
      _psource = s;
kpeter@605
   831
      _ptarget = t;
kpeter@605
   832
      _pstflow = k;
kpeter@605
   833
      return *this;
kpeter@605
   834
    }
kpeter@609
   835
    
kpeter@609
   836
    /// \brief Set the problem type.
kpeter@609
   837
    ///
kpeter@609
   838
    /// This function sets the problem type for the algorithm.
kpeter@609
   839
    /// If it is not used before calling \ref run(), the \ref GEQ problem
kpeter@609
   840
    /// type will be used.
kpeter@609
   841
    ///
kpeter@609
   842
    /// For more information see \ref ProblemType.
kpeter@609
   843
    ///
kpeter@609
   844
    /// \return <tt>(*this)</tt>
kpeter@609
   845
    NetworkSimplex& problemType(ProblemType problem_type) {
kpeter@609
   846
      _ptype = problem_type;
kpeter@609
   847
      return *this;
kpeter@609
   848
    }
kpeter@605
   849
kpeter@601
   850
    /// \brief Set the flow map.
kpeter@601
   851
    ///
kpeter@601
   852
    /// This function sets the flow map.
kpeter@605
   853
    /// If it is not used before calling \ref run(), an instance will
kpeter@605
   854
    /// be allocated automatically. The destructor deallocates this
kpeter@605
   855
    /// automatically allocated map, of course.
kpeter@601
   856
    ///
kpeter@601
   857
    /// \return <tt>(*this)</tt>
kpeter@605
   858
    NetworkSimplex& flowMap(FlowMap& map) {
kpeter@601
   859
      if (_local_flow) {
kpeter@603
   860
        delete _flow_map;
kpeter@601
   861
        _local_flow = false;
kpeter@601
   862
      }
kpeter@603
   863
      _flow_map = &map;
kpeter@601
   864
      return *this;
kpeter@601
   865
    }
kpeter@601
   866
kpeter@601
   867
    /// \brief Set the potential map.
kpeter@601
   868
    ///
kpeter@605
   869
    /// This function sets the potential map, which is used for storing
kpeter@605
   870
    /// the dual solution.
kpeter@605
   871
    /// If it is not used before calling \ref run(), an instance will
kpeter@605
   872
    /// be allocated automatically. The destructor deallocates this
kpeter@605
   873
    /// automatically allocated map, of course.
kpeter@601
   874
    ///
kpeter@601
   875
    /// \return <tt>(*this)</tt>
kpeter@605
   876
    NetworkSimplex& potentialMap(PotentialMap& map) {
kpeter@601
   877
      if (_local_potential) {
kpeter@603
   878
        delete _potential_map;
kpeter@601
   879
        _local_potential = false;
kpeter@601
   880
      }
kpeter@603
   881
      _potential_map = &map;
kpeter@601
   882
      return *this;
kpeter@601
   883
    }
kpeter@609
   884
    
kpeter@609
   885
    /// @}
kpeter@601
   886
kpeter@605
   887
    /// \name Execution Control
kpeter@605
   888
    /// The algorithm can be executed using \ref run().
kpeter@605
   889
kpeter@601
   890
    /// @{
kpeter@601
   891
kpeter@601
   892
    /// \brief Run the algorithm.
kpeter@601
   893
    ///
kpeter@601
   894
    /// This function runs the algorithm.
kpeter@609
   895
    /// The paramters can be specified using functions \ref lowerMap(),
kpeter@606
   896
    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
kpeter@609
   897
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
kpeter@609
   898
    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
kpeter@609
   899
    /// For example,
kpeter@605
   900
    /// \code
kpeter@605
   901
    ///   NetworkSimplex<ListDigraph> ns(graph);
kpeter@605
   902
    ///   ns.boundMaps(lower, upper).costMap(cost)
kpeter@605
   903
    ///     .supplyMap(sup).run();
kpeter@605
   904
    /// \endcode
kpeter@601
   905
    ///
kpeter@606
   906
    /// This function can be called more than once. All the parameters
kpeter@606
   907
    /// that have been given are kept for the next call, unless
kpeter@606
   908
    /// \ref reset() is called, thus only the modified parameters
kpeter@606
   909
    /// have to be set again. See \ref reset() for examples.
kpeter@606
   910
    ///
kpeter@605
   911
    /// \param pivot_rule The pivot rule that will be used during the
kpeter@605
   912
    /// algorithm. For more information see \ref PivotRule.
kpeter@601
   913
    ///
kpeter@601
   914
    /// \return \c true if a feasible flow can be found.
kpeter@605
   915
    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
kpeter@601
   916
      return init() && start(pivot_rule);
kpeter@601
   917
    }
kpeter@601
   918
kpeter@606
   919
    /// \brief Reset all the parameters that have been given before.
kpeter@606
   920
    ///
kpeter@606
   921
    /// This function resets all the paramaters that have been given
kpeter@609
   922
    /// before using functions \ref lowerMap(), \ref upperMap(),
kpeter@609
   923
    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
kpeter@609
   924
    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
kpeter@609
   925
    /// \ref flowMap() and \ref potentialMap().
kpeter@606
   926
    ///
kpeter@606
   927
    /// It is useful for multiple run() calls. If this function is not
kpeter@606
   928
    /// used, all the parameters given before are kept for the next
kpeter@606
   929
    /// \ref run() call.
kpeter@606
   930
    ///
kpeter@606
   931
    /// For example,
kpeter@606
   932
    /// \code
kpeter@606
   933
    ///   NetworkSimplex<ListDigraph> ns(graph);
kpeter@606
   934
    ///
kpeter@606
   935
    ///   // First run
kpeter@606
   936
    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
kpeter@606
   937
    ///     .supplyMap(sup).run();
kpeter@606
   938
    ///
kpeter@606
   939
    ///   // Run again with modified cost map (reset() is not called,
kpeter@606
   940
    ///   // so only the cost map have to be set again)
kpeter@606
   941
    ///   cost[e] += 100;
kpeter@606
   942
    ///   ns.costMap(cost).run();
kpeter@606
   943
    ///
kpeter@606
   944
    ///   // Run again from scratch using reset()
kpeter@606
   945
    ///   // (the lower bounds will be set to zero on all arcs)
kpeter@606
   946
    ///   ns.reset();
kpeter@606
   947
    ///   ns.capacityMap(cap).costMap(cost)
kpeter@606
   948
    ///     .supplyMap(sup).run();
kpeter@606
   949
    /// \endcode
kpeter@606
   950
    ///
kpeter@606
   951
    /// \return <tt>(*this)</tt>
kpeter@606
   952
    NetworkSimplex& reset() {
kpeter@606
   953
      delete _plower;
kpeter@606
   954
      delete _pupper;
kpeter@606
   955
      delete _pcost;
kpeter@606
   956
      delete _psupply;
kpeter@606
   957
      _plower = NULL;
kpeter@606
   958
      _pupper = NULL;
kpeter@606
   959
      _pcost = NULL;
kpeter@606
   960
      _psupply = NULL;
kpeter@606
   961
      _pstsup = false;
kpeter@609
   962
      _ptype = GEQ;
kpeter@609
   963
      if (_local_flow) delete _flow_map;
kpeter@609
   964
      if (_local_potential) delete _potential_map;
kpeter@609
   965
      _flow_map = NULL;
kpeter@609
   966
      _potential_map = NULL;
kpeter@609
   967
      _local_flow = false;
kpeter@609
   968
      _local_potential = false;
kpeter@609
   969
kpeter@606
   970
      return *this;
kpeter@606
   971
    }
kpeter@606
   972
kpeter@601
   973
    /// @}
kpeter@601
   974
kpeter@601
   975
    /// \name Query Functions
kpeter@601
   976
    /// The results of the algorithm can be obtained using these
kpeter@601
   977
    /// functions.\n
kpeter@605
   978
    /// The \ref run() function must be called before using them.
kpeter@605
   979
kpeter@601
   980
    /// @{
kpeter@601
   981
kpeter@605
   982
    /// \brief Return the total cost of the found flow.
kpeter@605
   983
    ///
kpeter@605
   984
    /// This function returns the total cost of the found flow.
kpeter@607
   985
    /// The complexity of the function is O(e).
kpeter@605
   986
    ///
kpeter@605
   987
    /// \note The return type of the function can be specified as a
kpeter@605
   988
    /// template parameter. For example,
kpeter@605
   989
    /// \code
kpeter@605
   990
    ///   ns.totalCost<double>();
kpeter@605
   991
    /// \endcode
kpeter@607
   992
    /// It is useful if the total cost cannot be stored in the \c Cost
kpeter@605
   993
    /// type of the algorithm, which is the default return type of the
kpeter@605
   994
    /// function.
kpeter@605
   995
    ///
kpeter@605
   996
    /// \pre \ref run() must be called before using this function.
kpeter@605
   997
    template <typename Num>
kpeter@605
   998
    Num totalCost() const {
kpeter@605
   999
      Num c = 0;
kpeter@605
  1000
      if (_pcost) {
kpeter@605
  1001
        for (ArcIt e(_graph); e != INVALID; ++e)
kpeter@605
  1002
          c += (*_flow_map)[e] * (*_pcost)[e];
kpeter@605
  1003
      } else {
kpeter@605
  1004
        for (ArcIt e(_graph); e != INVALID; ++e)
kpeter@605
  1005
          c += (*_flow_map)[e];
kpeter@605
  1006
      }
kpeter@605
  1007
      return c;
kpeter@605
  1008
    }
kpeter@605
  1009
kpeter@605
  1010
#ifndef DOXYGEN
kpeter@607
  1011
    Cost totalCost() const {
kpeter@607
  1012
      return totalCost<Cost>();
kpeter@605
  1013
    }
kpeter@605
  1014
#endif
kpeter@605
  1015
kpeter@605
  1016
    /// \brief Return the flow on the given arc.
kpeter@605
  1017
    ///
kpeter@605
  1018
    /// This function returns the flow on the given arc.
kpeter@605
  1019
    ///
kpeter@605
  1020
    /// \pre \ref run() must be called before using this function.
kpeter@607
  1021
    Flow flow(const Arc& a) const {
kpeter@605
  1022
      return (*_flow_map)[a];
kpeter@605
  1023
    }
kpeter@605
  1024
kpeter@601
  1025
    /// \brief Return a const reference to the flow map.
kpeter@601
  1026
    ///
kpeter@601
  1027
    /// This function returns a const reference to an arc map storing
kpeter@601
  1028
    /// the found flow.
kpeter@601
  1029
    ///
kpeter@601
  1030
    /// \pre \ref run() must be called before using this function.
kpeter@601
  1031
    const FlowMap& flowMap() const {
kpeter@603
  1032
      return *_flow_map;
kpeter@601
  1033
    }
kpeter@601
  1034
kpeter@605
  1035
    /// \brief Return the potential (dual value) of the given node.
kpeter@605
  1036
    ///
kpeter@605
  1037
    /// This function returns the potential (dual value) of the
kpeter@605
  1038
    /// given node.
kpeter@605
  1039
    ///
kpeter@605
  1040
    /// \pre \ref run() must be called before using this function.
kpeter@607
  1041
    Cost potential(const Node& n) const {
kpeter@605
  1042
      return (*_potential_map)[n];
kpeter@605
  1043
    }
kpeter@605
  1044
kpeter@601
  1045
    /// \brief Return a const reference to the potential map
kpeter@601
  1046
    /// (the dual solution).
kpeter@601
  1047
    ///
kpeter@601
  1048
    /// This function returns a const reference to a node map storing
kpeter@605
  1049
    /// the found potentials, which form the dual solution of the
kpeter@605
  1050
    /// \ref min_cost_flow "minimum cost flow" problem.
kpeter@601
  1051
    ///
kpeter@601
  1052
    /// \pre \ref run() must be called before using this function.
kpeter@601
  1053
    const PotentialMap& potentialMap() const {
kpeter@603
  1054
      return *_potential_map;
kpeter@601
  1055
    }
kpeter@601
  1056
kpeter@601
  1057
    /// @}
kpeter@601
  1058
kpeter@601
  1059
  private:
kpeter@601
  1060
kpeter@601
  1061
    // Initialize internal data structures
kpeter@601
  1062
    bool init() {
kpeter@601
  1063
      // Initialize result maps
kpeter@603
  1064
      if (!_flow_map) {
kpeter@603
  1065
        _flow_map = new FlowMap(_graph);
kpeter@601
  1066
        _local_flow = true;
kpeter@601
  1067
      }
kpeter@603
  1068
      if (!_potential_map) {
kpeter@603
  1069
        _potential_map = new PotentialMap(_graph);
kpeter@601
  1070
        _local_potential = true;
kpeter@601
  1071
      }
kpeter@601
  1072
kpeter@601
  1073
      // Initialize vectors
kpeter@603
  1074
      _node_num = countNodes(_graph);
kpeter@603
  1075
      _arc_num = countArcs(_graph);
kpeter@601
  1076
      int all_node_num = _node_num + 1;
kpeter@603
  1077
      int all_arc_num = _arc_num + _node_num;
kpeter@605
  1078
      if (_node_num == 0) return false;
kpeter@601
  1079
kpeter@603
  1080
      _arc_ref.resize(_arc_num);
kpeter@603
  1081
      _source.resize(all_arc_num);
kpeter@603
  1082
      _target.resize(all_arc_num);
kpeter@601
  1083
kpeter@603
  1084
      _cap.resize(all_arc_num);
kpeter@603
  1085
      _cost.resize(all_arc_num);
kpeter@601
  1086
      _supply.resize(all_node_num);
kpeter@606
  1087
      _flow.resize(all_arc_num);
kpeter@606
  1088
      _pi.resize(all_node_num);
kpeter@601
  1089
kpeter@601
  1090
      _parent.resize(all_node_num);
kpeter@601
  1091
      _pred.resize(all_node_num);
kpeter@603
  1092
      _forward.resize(all_node_num);
kpeter@601
  1093
      _thread.resize(all_node_num);
kpeter@604
  1094
      _rev_thread.resize(all_node_num);
kpeter@604
  1095
      _succ_num.resize(all_node_num);
kpeter@604
  1096
      _last_succ.resize(all_node_num);
kpeter@606
  1097
      _state.resize(all_arc_num);
kpeter@601
  1098
kpeter@601
  1099
      // Initialize node related data
kpeter@601
  1100
      bool valid_supply = true;
kpeter@609
  1101
      Flow sum_supply = 0;
kpeter@605
  1102
      if (!_pstsup && !_psupply) {
kpeter@605
  1103
        _pstsup = true;
kpeter@605
  1104
        _psource = _ptarget = NodeIt(_graph);
kpeter@605
  1105
        _pstflow = 0;
kpeter@605
  1106
      }
kpeter@605
  1107
      if (_psupply) {
kpeter@601
  1108
        int i = 0;
kpeter@603
  1109
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
kpeter@601
  1110
          _node_id[n] = i;
kpeter@605
  1111
          _supply[i] = (*_psupply)[n];
kpeter@609
  1112
          sum_supply += _supply[i];
kpeter@601
  1113
        }
kpeter@609
  1114
        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
kpeter@609
  1115
                       (_ptype == LEQ && sum_supply >= 0);
kpeter@601
  1116
      } else {
kpeter@601
  1117
        int i = 0;
kpeter@603
  1118
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
kpeter@601
  1119
          _node_id[n] = i;
kpeter@601
  1120
          _supply[i] = 0;
kpeter@601
  1121
        }
kpeter@605
  1122
        _supply[_node_id[_psource]] =  _pstflow;
kpeter@609
  1123
        _supply[_node_id[_ptarget]] = -_pstflow;
kpeter@601
  1124
      }
kpeter@601
  1125
      if (!valid_supply) return false;
kpeter@601
  1126
kpeter@609
  1127
      // Infinite capacity value
kpeter@609
  1128
      Flow inf_cap =
kpeter@609
  1129
        std::numeric_limits<Flow>::has_infinity ?
kpeter@609
  1130
        std::numeric_limits<Flow>::infinity() :
kpeter@609
  1131
        std::numeric_limits<Flow>::max();
kpeter@609
  1132
kpeter@609
  1133
      // Initialize artifical cost
kpeter@609
  1134
      Cost art_cost;
kpeter@609
  1135
      if (std::numeric_limits<Cost>::is_exact) {
kpeter@609
  1136
        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
kpeter@609
  1137
      } else {
kpeter@609
  1138
        art_cost = std::numeric_limits<Cost>::min();
kpeter@609
  1139
        for (int i = 0; i != _arc_num; ++i) {
kpeter@609
  1140
          if (_cost[i] > art_cost) art_cost = _cost[i];
kpeter@609
  1141
        }
kpeter@609
  1142
        art_cost = (art_cost + 1) * _node_num;
kpeter@609
  1143
      }
kpeter@609
  1144
kpeter@609
  1145
      // Run Circulation to check if a feasible solution exists
kpeter@609
  1146
      typedef ConstMap<Arc, Flow> ConstArcMap;
kpeter@609
  1147
      FlowNodeMap *csup = NULL;
kpeter@609
  1148
      bool local_csup = false;
kpeter@609
  1149
      if (_psupply) {
kpeter@609
  1150
        csup = _psupply;
kpeter@609
  1151
      } else {
kpeter@609
  1152
        csup = new FlowNodeMap(_graph, 0);
kpeter@609
  1153
        (*csup)[_psource] =  _pstflow;
kpeter@609
  1154
        (*csup)[_ptarget] = -_pstflow;
kpeter@609
  1155
        local_csup = true;
kpeter@609
  1156
      }
kpeter@609
  1157
      bool circ_result = false;
kpeter@609
  1158
      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
kpeter@609
  1159
        // GEQ problem type
kpeter@609
  1160
        if (_plower) {
kpeter@609
  1161
          if (_pupper) {
kpeter@609
  1162
            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
kpeter@609
  1163
              circ(_graph, *_plower, *_pupper, *csup);
kpeter@609
  1164
            circ_result = circ.run();
kpeter@609
  1165
          } else {
kpeter@609
  1166
            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
kpeter@609
  1167
              circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
kpeter@609
  1168
            circ_result = circ.run();
kpeter@609
  1169
          }
kpeter@609
  1170
        } else {
kpeter@609
  1171
          if (_pupper) {
kpeter@609
  1172
            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
kpeter@609
  1173
              circ(_graph, ConstArcMap(0), *_pupper, *csup);
kpeter@609
  1174
            circ_result = circ.run();
kpeter@609
  1175
          } else {
kpeter@609
  1176
            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
kpeter@609
  1177
              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
kpeter@609
  1178
            circ_result = circ.run();
kpeter@609
  1179
          }
kpeter@609
  1180
        }
kpeter@609
  1181
      } else {
kpeter@609
  1182
        // LEQ problem type
kpeter@609
  1183
        typedef ReverseDigraph<const GR> RevGraph;
kpeter@609
  1184
        typedef NegMap<FlowNodeMap> NegNodeMap;
kpeter@609
  1185
        RevGraph rgraph(_graph);
kpeter@609
  1186
        NegNodeMap neg_csup(*csup);
kpeter@609
  1187
        if (_plower) {
kpeter@609
  1188
          if (_pupper) {
kpeter@609
  1189
            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
kpeter@609
  1190
              circ(rgraph, *_plower, *_pupper, neg_csup);
kpeter@609
  1191
            circ_result = circ.run();
kpeter@609
  1192
          } else {
kpeter@609
  1193
            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
kpeter@609
  1194
              circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
kpeter@609
  1195
            circ_result = circ.run();
kpeter@609
  1196
          }
kpeter@609
  1197
        } else {
kpeter@609
  1198
          if (_pupper) {
kpeter@609
  1199
            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
kpeter@609
  1200
              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
kpeter@609
  1201
            circ_result = circ.run();
kpeter@609
  1202
          } else {
kpeter@609
  1203
            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
kpeter@609
  1204
              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
kpeter@609
  1205
            circ_result = circ.run();
kpeter@609
  1206
          }
kpeter@609
  1207
        }
kpeter@609
  1208
      }
kpeter@609
  1209
      if (local_csup) delete csup;
kpeter@609
  1210
      if (!circ_result) return false;
kpeter@609
  1211
kpeter@601
  1212
      // Set data for the artificial root node
kpeter@601
  1213
      _root = _node_num;
kpeter@601
  1214
      _parent[_root] = -1;
kpeter@601
  1215
      _pred[_root] = -1;
kpeter@601
  1216
      _thread[_root] = 0;
kpeter@604
  1217
      _rev_thread[0] = _root;
kpeter@604
  1218
      _succ_num[_root] = all_node_num;
kpeter@604
  1219
      _last_succ[_root] = _root - 1;
kpeter@609
  1220
      _supply[_root] = -sum_supply;
kpeter@609
  1221
      if (sum_supply < 0) {
kpeter@609
  1222
        _pi[_root] = -art_cost;
kpeter@609
  1223
      } else {
kpeter@609
  1224
        _pi[_root] = art_cost;
kpeter@609
  1225
      }
kpeter@601
  1226
kpeter@601
  1227
      // Store the arcs in a mixed order
kpeter@601
  1228
      int k = std::max(int(sqrt(_arc_num)), 10);
kpeter@601
  1229
      int i = 0;
kpeter@603
  1230
      for (ArcIt e(_graph); e != INVALID; ++e) {
kpeter@603
  1231
        _arc_ref[i] = e;
kpeter@601
  1232
        if ((i += k) >= _arc_num) i = (i % k) + 1;
kpeter@601
  1233
      }
kpeter@601
  1234
kpeter@601
  1235
      // Initialize arc maps
kpeter@605
  1236
      if (_pupper && _pcost) {
kpeter@605
  1237
        for (int i = 0; i != _arc_num; ++i) {
kpeter@605
  1238
          Arc e = _arc_ref[i];
kpeter@605
  1239
          _source[i] = _node_id[_graph.source(e)];
kpeter@605
  1240
          _target[i] = _node_id[_graph.target(e)];
kpeter@605
  1241
          _cap[i] = (*_pupper)[e];
kpeter@605
  1242
          _cost[i] = (*_pcost)[e];
kpeter@606
  1243
          _flow[i] = 0;
kpeter@606
  1244
          _state[i] = STATE_LOWER;
kpeter@605
  1245
        }
kpeter@605
  1246
      } else {
kpeter@605
  1247
        for (int i = 0; i != _arc_num; ++i) {
kpeter@605
  1248
          Arc e = _arc_ref[i];
kpeter@605
  1249
          _source[i] = _node_id[_graph.source(e)];
kpeter@605
  1250
          _target[i] = _node_id[_graph.target(e)];
kpeter@606
  1251
          _flow[i] = 0;
kpeter@606
  1252
          _state[i] = STATE_LOWER;
kpeter@605
  1253
        }
kpeter@605
  1254
        if (_pupper) {
kpeter@605
  1255
          for (int i = 0; i != _arc_num; ++i)
kpeter@605
  1256
            _cap[i] = (*_pupper)[_arc_ref[i]];
kpeter@605
  1257
        } else {
kpeter@605
  1258
          for (int i = 0; i != _arc_num; ++i)
kpeter@608
  1259
            _cap[i] = inf_cap;
kpeter@605
  1260
        }
kpeter@605
  1261
        if (_pcost) {
kpeter@605
  1262
          for (int i = 0; i != _arc_num; ++i)
kpeter@605
  1263
            _cost[i] = (*_pcost)[_arc_ref[i]];
kpeter@605
  1264
        } else {
kpeter@605
  1265
          for (int i = 0; i != _arc_num; ++i)
kpeter@605
  1266
            _cost[i] = 1;
kpeter@605
  1267
        }
kpeter@601
  1268
      }
kpeter@608
  1269
      
kpeter@601
  1270
      // Remove non-zero lower bounds
kpeter@605
  1271
      if (_plower) {
kpeter@601
  1272
        for (int i = 0; i != _arc_num; ++i) {
kpeter@607
  1273
          Flow c = (*_plower)[_arc_ref[i]];
kpeter@601
  1274
          if (c != 0) {
kpeter@601
  1275
            _cap[i] -= c;
kpeter@601
  1276
            _supply[_source[i]] -= c;
kpeter@601
  1277
            _supply[_target[i]] += c;
kpeter@601
  1278
          }
kpeter@601
  1279
        }
kpeter@601
  1280
      }
kpeter@601
  1281
kpeter@601
  1282
      // Add artificial arcs and initialize the spanning tree data structure
kpeter@601
  1283
      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
kpeter@601
  1284
        _thread[u] = u + 1;
kpeter@604
  1285
        _rev_thread[u + 1] = u;
kpeter@604
  1286
        _succ_num[u] = 1;
kpeter@604
  1287
        _last_succ[u] = u;
kpeter@601
  1288
        _parent[u] = _root;
kpeter@601
  1289
        _pred[u] = e;
kpeter@608
  1290
        _cost[e] = art_cost;
kpeter@608
  1291
        _cap[e] = inf_cap;
kpeter@606
  1292
        _state[e] = STATE_TREE;
kpeter@609
  1293
        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
kpeter@601
  1294
          _flow[e] = _supply[u];
kpeter@601
  1295
          _forward[u] = true;
kpeter@609
  1296
          _pi[u] = -art_cost + _pi[_root];
kpeter@601
  1297
        } else {
kpeter@601
  1298
          _flow[e] = -_supply[u];
kpeter@601
  1299
          _forward[u] = false;
kpeter@609
  1300
          _pi[u] = art_cost + _pi[_root];
kpeter@601
  1301
        }
kpeter@601
  1302
      }
kpeter@601
  1303
kpeter@601
  1304
      return true;
kpeter@601
  1305
    }
kpeter@601
  1306
kpeter@601
  1307
    // Find the join node
kpeter@601
  1308
    void findJoinNode() {
kpeter@603
  1309
      int u = _source[in_arc];
kpeter@603
  1310
      int v = _target[in_arc];
kpeter@601
  1311
      while (u != v) {
kpeter@604
  1312
        if (_succ_num[u] < _succ_num[v]) {
kpeter@604
  1313
          u = _parent[u];
kpeter@604
  1314
        } else {
kpeter@604
  1315
          v = _parent[v];
kpeter@604
  1316
        }
kpeter@601
  1317
      }
kpeter@601
  1318
      join = u;
kpeter@601
  1319
    }
kpeter@601
  1320
kpeter@601
  1321
    // Find the leaving arc of the cycle and returns true if the
kpeter@601
  1322
    // leaving arc is not the same as the entering arc
kpeter@601
  1323
    bool findLeavingArc() {
kpeter@601
  1324
      // Initialize first and second nodes according to the direction
kpeter@601
  1325
      // of the cycle
kpeter@603
  1326
      if (_state[in_arc] == STATE_LOWER) {
kpeter@603
  1327
        first  = _source[in_arc];
kpeter@603
  1328
        second = _target[in_arc];
kpeter@601
  1329
      } else {
kpeter@603
  1330
        first  = _target[in_arc];
kpeter@603
  1331
        second = _source[in_arc];
kpeter@601
  1332
      }
kpeter@603
  1333
      delta = _cap[in_arc];
kpeter@601
  1334
      int result = 0;
kpeter@607
  1335
      Flow d;
kpeter@601
  1336
      int e;
kpeter@601
  1337
kpeter@601
  1338
      // Search the cycle along the path form the first node to the root
kpeter@601
  1339
      for (int u = first; u != join; u = _parent[u]) {
kpeter@601
  1340
        e = _pred[u];
kpeter@601
  1341
        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
kpeter@601
  1342
        if (d < delta) {
kpeter@601
  1343
          delta = d;
kpeter@601
  1344
          u_out = u;
kpeter@601
  1345
          result = 1;
kpeter@601
  1346
        }
kpeter@601
  1347
      }
kpeter@601
  1348
      // Search the cycle along the path form the second node to the root
kpeter@601
  1349
      for (int u = second; u != join; u = _parent[u]) {
kpeter@601
  1350
        e = _pred[u];
kpeter@601
  1351
        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
kpeter@601
  1352
        if (d <= delta) {
kpeter@601
  1353
          delta = d;
kpeter@601
  1354
          u_out = u;
kpeter@601
  1355
          result = 2;
kpeter@601
  1356
        }
kpeter@601
  1357
      }
kpeter@601
  1358
kpeter@601
  1359
      if (result == 1) {
kpeter@601
  1360
        u_in = first;
kpeter@601
  1361
        v_in = second;
kpeter@601
  1362
      } else {
kpeter@601
  1363
        u_in = second;
kpeter@601
  1364
        v_in = first;
kpeter@601
  1365
      }
kpeter@601
  1366
      return result != 0;
kpeter@601
  1367
    }
kpeter@601
  1368
kpeter@601
  1369
    // Change _flow and _state vectors
kpeter@601
  1370
    void changeFlow(bool change) {
kpeter@601
  1371
      // Augment along the cycle
kpeter@601
  1372
      if (delta > 0) {
kpeter@607
  1373
        Flow val = _state[in_arc] * delta;
kpeter@603
  1374
        _flow[in_arc] += val;
kpeter@603
  1375
        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
kpeter@601
  1376
          _flow[_pred[u]] += _forward[u] ? -val : val;
kpeter@601
  1377
        }
kpeter@603
  1378
        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
kpeter@601
  1379
          _flow[_pred[u]] += _forward[u] ? val : -val;
kpeter@601
  1380
        }
kpeter@601
  1381
      }
kpeter@601
  1382
      // Update the state of the entering and leaving arcs
kpeter@601
  1383
      if (change) {
kpeter@603
  1384
        _state[in_arc] = STATE_TREE;
kpeter@601
  1385
        _state[_pred[u_out]] =
kpeter@601
  1386
          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
kpeter@601
  1387
      } else {
kpeter@603
  1388
        _state[in_arc] = -_state[in_arc];
kpeter@601
  1389
      }
kpeter@601
  1390
    }
kpeter@601
  1391
kpeter@604
  1392
    // Update the tree structure
kpeter@604
  1393
    void updateTreeStructure() {
kpeter@604
  1394
      int u, w;
kpeter@604
  1395
      int old_rev_thread = _rev_thread[u_out];
kpeter@604
  1396
      int old_succ_num = _succ_num[u_out];
kpeter@604
  1397
      int old_last_succ = _last_succ[u_out];
kpeter@601
  1398
      v_out = _parent[u_out];
kpeter@601
  1399
kpeter@604
  1400
      u = _last_succ[u_in];  // the last successor of u_in
kpeter@604
  1401
      right = _thread[u];    // the node after it
kpeter@604
  1402
kpeter@604
  1403
      // Handle the case when old_rev_thread equals to v_in
kpeter@604
  1404
      // (it also means that join and v_out coincide)
kpeter@604
  1405
      if (old_rev_thread == v_in) {
kpeter@604
  1406
        last = _thread[_last_succ[u_out]];
kpeter@604
  1407
      } else {
kpeter@604
  1408
        last = _thread[v_in];
kpeter@601
  1409
      }
kpeter@601
  1410
kpeter@604
  1411
      // Update _thread and _parent along the stem nodes (i.e. the nodes
kpeter@604
  1412
      // between u_in and u_out, whose parent have to be changed)
kpeter@601
  1413
      _thread[v_in] = stem = u_in;
kpeter@604
  1414
      _dirty_revs.clear();
kpeter@604
  1415
      _dirty_revs.push_back(v_in);
kpeter@601
  1416
      par_stem = v_in;
kpeter@601
  1417
      while (stem != u_out) {
kpeter@604
  1418
        // Insert the next stem node into the thread list
kpeter@604
  1419
        new_stem = _parent[stem];
kpeter@604
  1420
        _thread[u] = new_stem;
kpeter@604
  1421
        _dirty_revs.push_back(u);
kpeter@601
  1422
kpeter@604
  1423
        // Remove the subtree of stem from the thread list
kpeter@604
  1424
        w = _rev_thread[stem];
kpeter@604
  1425
        _thread[w] = right;
kpeter@604
  1426
        _rev_thread[right] = w;
kpeter@601
  1427
kpeter@604
  1428
        // Change the parent node and shift stem nodes
kpeter@601
  1429
        _parent[stem] = par_stem;
kpeter@601
  1430
        par_stem = stem;
kpeter@601
  1431
        stem = new_stem;
kpeter@601
  1432
kpeter@604
  1433
        // Update u and right
kpeter@604
  1434
        u = _last_succ[stem] == _last_succ[par_stem] ?
kpeter@604
  1435
          _rev_thread[par_stem] : _last_succ[stem];
kpeter@601
  1436
        right = _thread[u];
kpeter@601
  1437
      }
kpeter@601
  1438
      _parent[u_out] = par_stem;
kpeter@601
  1439
      _thread[u] = last;
kpeter@604
  1440
      _rev_thread[last] = u;
kpeter@604
  1441
      _last_succ[u_out] = u;
kpeter@601
  1442
kpeter@604
  1443
      // Remove the subtree of u_out from the thread list except for
kpeter@604
  1444
      // the case when old_rev_thread equals to v_in
kpeter@604
  1445
      // (it also means that join and v_out coincide)
kpeter@604
  1446
      if (old_rev_thread != v_in) {
kpeter@604
  1447
        _thread[old_rev_thread] = right;
kpeter@604
  1448
        _rev_thread[right] = old_rev_thread;
kpeter@604
  1449
      }
kpeter@604
  1450
kpeter@604
  1451
      // Update _rev_thread using the new _thread values
kpeter@604
  1452
      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
kpeter@604
  1453
        u = _dirty_revs[i];
kpeter@604
  1454
        _rev_thread[_thread[u]] = u;
kpeter@604
  1455
      }
kpeter@604
  1456
kpeter@604
  1457
      // Update _pred, _forward, _last_succ and _succ_num for the
kpeter@604
  1458
      // stem nodes from u_out to u_in
kpeter@604
  1459
      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
kpeter@604
  1460
      u = u_out;
kpeter@604
  1461
      while (u != u_in) {
kpeter@604
  1462
        w = _parent[u];
kpeter@604
  1463
        _pred[u] = _pred[w];
kpeter@604
  1464
        _forward[u] = !_forward[w];
kpeter@604
  1465
        tmp_sc += _succ_num[u] - _succ_num[w];
kpeter@604
  1466
        _succ_num[u] = tmp_sc;
kpeter@604
  1467
        _last_succ[w] = tmp_ls;
kpeter@604
  1468
        u = w;
kpeter@604
  1469
      }
kpeter@604
  1470
      _pred[u_in] = in_arc;
kpeter@604
  1471
      _forward[u_in] = (u_in == _source[in_arc]);
kpeter@604
  1472
      _succ_num[u_in] = old_succ_num;
kpeter@604
  1473
kpeter@604
  1474
      // Set limits for updating _last_succ form v_in and v_out
kpeter@604
  1475
      // towards the root
kpeter@604
  1476
      int up_limit_in = -1;
kpeter@604
  1477
      int up_limit_out = -1;
kpeter@604
  1478
      if (_last_succ[join] == v_in) {
kpeter@604
  1479
        up_limit_out = join;
kpeter@601
  1480
      } else {
kpeter@604
  1481
        up_limit_in = join;
kpeter@604
  1482
      }
kpeter@604
  1483
kpeter@604
  1484
      // Update _last_succ from v_in towards the root
kpeter@604
  1485
      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
kpeter@604
  1486
           u = _parent[u]) {
kpeter@604
  1487
        _last_succ[u] = _last_succ[u_out];
kpeter@604
  1488
      }
kpeter@604
  1489
      // Update _last_succ from v_out towards the root
kpeter@604
  1490
      if (join != old_rev_thread && v_in != old_rev_thread) {
kpeter@604
  1491
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
kpeter@604
  1492
             u = _parent[u]) {
kpeter@604
  1493
          _last_succ[u] = old_rev_thread;
kpeter@604
  1494
        }
kpeter@604
  1495
      } else {
kpeter@604
  1496
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
kpeter@604
  1497
             u = _parent[u]) {
kpeter@604
  1498
          _last_succ[u] = _last_succ[u_out];
kpeter@604
  1499
        }
kpeter@604
  1500
      }
kpeter@604
  1501
kpeter@604
  1502
      // Update _succ_num from v_in to join
kpeter@604
  1503
      for (u = v_in; u != join; u = _parent[u]) {
kpeter@604
  1504
        _succ_num[u] += old_succ_num;
kpeter@604
  1505
      }
kpeter@604
  1506
      // Update _succ_num from v_out to join
kpeter@604
  1507
      for (u = v_out; u != join; u = _parent[u]) {
kpeter@604
  1508
        _succ_num[u] -= old_succ_num;
kpeter@601
  1509
      }
kpeter@601
  1510
    }
kpeter@601
  1511
kpeter@604
  1512
    // Update potentials
kpeter@604
  1513
    void updatePotential() {
kpeter@607
  1514
      Cost sigma = _forward[u_in] ?
kpeter@601
  1515
        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
kpeter@601
  1516
        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
kpeter@608
  1517
      // Update potentials in the subtree, which has been moved
kpeter@608
  1518
      int end = _thread[_last_succ[u_in]];
kpeter@608
  1519
      for (int u = u_in; u != end; u = _thread[u]) {
kpeter@608
  1520
        _pi[u] += sigma;
kpeter@601
  1521
      }
kpeter@601
  1522
    }
kpeter@601
  1523
kpeter@601
  1524
    // Execute the algorithm
kpeter@605
  1525
    bool start(PivotRule pivot_rule) {
kpeter@601
  1526
      // Select the pivot rule implementation
kpeter@601
  1527
      switch (pivot_rule) {
kpeter@605
  1528
        case FIRST_ELIGIBLE:
kpeter@601
  1529
          return start<FirstEligiblePivotRule>();
kpeter@605
  1530
        case BEST_ELIGIBLE:
kpeter@601
  1531
          return start<BestEligiblePivotRule>();
kpeter@605
  1532
        case BLOCK_SEARCH:
kpeter@601
  1533
          return start<BlockSearchPivotRule>();
kpeter@605
  1534
        case CANDIDATE_LIST:
kpeter@601
  1535
          return start<CandidateListPivotRule>();
kpeter@605
  1536
        case ALTERING_LIST:
kpeter@601
  1537
          return start<AlteringListPivotRule>();
kpeter@601
  1538
      }
kpeter@601
  1539
      return false;
kpeter@601
  1540
    }
kpeter@601
  1541
kpeter@605
  1542
    template <typename PivotRuleImpl>
kpeter@601
  1543
    bool start() {
kpeter@605
  1544
      PivotRuleImpl pivot(*this);
kpeter@601
  1545
kpeter@605
  1546
      // Execute the Network Simplex algorithm
kpeter@601
  1547
      while (pivot.findEnteringArc()) {
kpeter@601
  1548
        findJoinNode();
kpeter@601
  1549
        bool change = findLeavingArc();
kpeter@601
  1550
        changeFlow(change);
kpeter@601
  1551
        if (change) {
kpeter@604
  1552
          updateTreeStructure();
kpeter@604
  1553
          updatePotential();
kpeter@601
  1554
        }
kpeter@601
  1555
      }
kpeter@601
  1556
kpeter@603
  1557
      // Copy flow values to _flow_map
kpeter@605
  1558
      if (_plower) {
kpeter@601
  1559
        for (int i = 0; i != _arc_num; ++i) {
kpeter@603
  1560
          Arc e = _arc_ref[i];
kpeter@605
  1561
          _flow_map->set(e, (*_plower)[e] + _flow[i]);
kpeter@601
  1562
        }
kpeter@601
  1563
      } else {
kpeter@601
  1564
        for (int i = 0; i != _arc_num; ++i) {
kpeter@603
  1565
          _flow_map->set(_arc_ref[i], _flow[i]);
kpeter@601
  1566
        }
kpeter@601
  1567
      }
kpeter@603
  1568
      // Copy potential values to _potential_map
kpeter@603
  1569
      for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@603
  1570
        _potential_map->set(n, _pi[_node_id[n]]);
kpeter@601
  1571
      }
kpeter@601
  1572
kpeter@601
  1573
      return true;
kpeter@601
  1574
    }
kpeter@601
  1575
kpeter@601
  1576
  }; //class NetworkSimplex
kpeter@601
  1577
kpeter@601
  1578
  ///@}
kpeter@601
  1579
kpeter@601
  1580
} //namespace lemon
kpeter@601
  1581
kpeter@601
  1582
#endif //LEMON_NETWORK_SIMPLEX_H