RhodeCode

EXTRA_DIST += \

	doc/Doxyfile.in \

	doc/DoxygenLayout.xml \

	doc/coding_style.dox \

	doc/dirs.dox \

	doc/groups.dox \

	doc/lgf.dox \

	doc/license.dox \

	doc/mainpage.dox \

	doc/migration.dox \

	doc/named-param.dox \

	doc/namespaces.dox \

	doc/html \

	doc/CMakeLists.txt

DOC_EPS_IMAGES18 = \

	grid_graph.eps \

	nodeshape_0.eps \

	nodeshape_1.eps \

	nodeshape_2.eps \

	nodeshape_3.eps \

	nodeshape_4.eps

DOC_EPS_IMAGES27 = \

	bipartite_matching.eps \

	bipartite_partitions.eps \

- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm.

- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees.

- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees.

In most cases the \ref Preflow "Preflow" algorithm provides the

fastest method for computing a maximum flow. All implementations

also provide functions to query the minimum cut, which is the dual

problem of maximum flow.

\ref Circulation is a preflow push-relabel algorithm implemented directly 

for finding feasible circulations, which is a somewhat different problem,

but it is strongly related to maximum flow.

For more information, see \ref Circulation.

*/

/**

@ingroup algs

\brief Algorithms for finding minimum cost flows and circulations.

This group contains the algorithms for finding minimum cost flows and

 - \ref NetworkSimplex Primal Network Simplex algorithm with various

   pivot strategies.

 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on

   cost scaling.

 - \ref CapacityScaling Successive Shortest %Path algorithm with optional

   capacity scaling.

 - \ref CancelAndTighten The Cancel and Tighten algorithm.

 - \ref CycleCanceling Cycle-Canceling algorithms.

In general NetworkSimplex is the most efficient implementation,

but in special cases other algorithms could be faster.

For example, if the total supply and/or capacities are rather small,

CapacityScaling is usually the fastest algorithm (without effective scaling).

*/

/**

@defgroup min_cut Minimum Cut Algorithms

@ingroup algs

\brief Algorithms for finding minimum cut in graphs.

This group contains the algorithms for finding minimum cut in graphs.

The \e minimum \e cut \e problem is to find a non-empty and non-complete

\f$X\f$ subset of the nodes with minimum overall capacity on

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_NETWORK_SIMPLEX_H

#define LEMON_NETWORK_SIMPLEX_H

///

/// \file

/// \brief Network Simplex algorithm for finding a minimum cost flow.

#include <vector>

#include <limits>

#include <algorithm>

#include <lemon/core.h>

#include <lemon/math.h>

namespace lemon {

  /// @{

  /// \brief Implementation of the primal Network Simplex algorithm

  /// for finding a \ref min_cost_flow "minimum cost flow".

  ///

  /// \ref NetworkSimplex implements the primal Network Simplex algorithm

  /// for finding a \ref min_cost_flow "minimum cost flow".

  /// This algorithm is a specialized version of the linear programming

  /// simplex method directly for the minimum cost flow problem.

  /// It is one of the most efficient solution methods.

  ///

  /// In general this class is the fastest implementation available

  /// in LEMON for the minimum cost flow problem.

  /// Moreover it supports both directions of the supply/demand inequality

  /// constraints. For more information see \ref SupplyType.

  ///

      /// The problem has optimal solution (i.e. it is feasible and

      /// bounded), and the algorithm has found optimal flow and node

      /// potentials (primal and dual solutions).

      OPTIMAL,

      /// The objective function of the problem is unbounded, i.e.

      /// there is a directed cycle having negative total cost and

      /// infinite upper bound.

      UNBOUNDED

    };

    /// \brief Constants for selecting the type of the supply constraints.

    ///

    /// Enum type containing constants for selecting the supply type,

    /// i.e. the direction of the inequalities in the supply/demand

    /// constraints of the \ref min_cost_flow "minimum cost flow problem".

    ///

    enum SupplyType {

      /// This option means that there are <em>"greater or equal"</em>

      GEQ,

      /// This option means that there are <em>"less or equal"</em>

    };

    /// \brief Constants for selecting the pivot rule.

    ///

    /// Enum type containing constants for selecting the pivot rule for

    /// the \ref run() function.

    ///

    /// \ref NetworkSimplex provides five different pivot rule

    /// implementations that significantly affect the running time

    /// of the algorithm.

    /// By default \ref BLOCK_SEARCH "Block Search" is used, which

    /// proved to be the most efficient and the most robust on various

    /// test inputs according to our benchmark tests.

    /// However another pivot rule can be selected using the \ref run()

    /// function with the proper parameter.

    enum PivotRule {

    typedef std::vector<Value> ValueVector;

    typedef std::vector<Cost> CostVector;

    // State constants for arcs

    enum ArcStateEnum {

      STATE_UPPER = -1,

      STATE_TREE  =  0,

      STATE_LOWER =  1

    };

  private:

    // Data related to the underlying digraph

    const GR &_graph;

    int _node_num;

    int _arc_num;

    // Parameters of the problem

    bool _have_lower;

    SupplyType _stype;

    Value _sum_supply;

    // Data structures for storing the digraph

    IntNodeMap _node_id;

    IntArcMap _arc_id;

    IntVector _source;

    IntVector _target;

    // Node and arc data

    ValueVector _lower;

    ValueVector _upper;

    ValueVector _cap;

    const Value INF;

  private:

    // Implementation of the First Eligible pivot rule

    class FirstEligiblePivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const IntVector  &_state;

      const CostVector &_pi;

      int &_in_arc;

      // Pivot rule data

      int _next_arc;

    public:

      // Constructor

      FirstEligiblePivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

      {}

      // Find next entering arc

      bool findEnteringArc() {

        Cost c;

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < 0) {

            _in_arc = e;

            _next_arc = e + 1;

            return true;

          }

        }

        for (int e = 0; e < _next_arc; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < 0) {

            _in_arc = e;

            _next_arc = e + 1;

            return true;

          }

        }

        return false;

    }; //class FirstEligiblePivotRule

    // Implementation of the Best Eligible pivot rule

    class BestEligiblePivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const IntVector  &_state;

      const CostVector &_pi;

      int &_in_arc;

    public:

      // Constructor

      BestEligiblePivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

      {}

      // Find next entering arc

      bool findEnteringArc() {

        Cost c, min = 0;

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < min) {

            min = c;

            _in_arc = e;

          }

        }

        return min < 0;

      }

    }; //class BestEligiblePivotRule

    // Implementation of the Block Search pivot rule

    class BlockSearchPivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const IntVector  &_state;

      const CostVector &_pi;

      int &_in_arc;

      // Pivot rule data

      int _block_size;

      int _next_arc;

    public:

      // Constructor

      BlockSearchPivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

      {

        // The main parameters of the pivot rule

        const int MIN_BLOCK_SIZE = 10;

        _block_size = std::max( int(BLOCK_SIZE_FACTOR *

                                MIN_BLOCK_SIZE );

      }

      // Find next entering arc

      bool findEnteringArc() {

        Cost c, min = 0;

        int cnt = _block_size;

        int e, min_arc = _next_arc;

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < min) {

            min = c;

            min_arc = e;

          }

          if (--cnt == 0) {

            if (min < 0) break;

            cnt = _block_size;

          }

        }

        if (min == 0 || cnt > 0) {

          for (e = 0; e < _next_arc; ++e) {

            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

            if (c < min) {

              min = c;

              min_arc = e;

    }; //class BlockSearchPivotRule

    // Implementation of the Candidate List pivot rule

    class CandidateListPivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const IntVector  &_state;

      const CostVector &_pi;

      int &_in_arc;

      // Pivot rule data

      IntVector _candidates;

      int _list_length, _minor_limit;

      int _curr_length, _minor_count;

      int _next_arc;

    public:

      /// Constructor

      CandidateListPivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

      {

        // The main parameters of the pivot rule

        const double LIST_LENGTH_FACTOR = 1.0;

        const int MIN_LIST_LENGTH = 10;

        const double MINOR_LIMIT_FACTOR = 0.1;

        const int MIN_MINOR_LIMIT = 3;

        _list_length = std::max( int(LIST_LENGTH_FACTOR *

                                 MIN_LIST_LENGTH );

        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),

                                 MIN_MINOR_LIMIT );

        _curr_length = _minor_count = 0;

        _candidates.resize(_list_length);

      }

      /// Find next entering arc

      bool findEnteringArc() {

        Cost min, c;

        int e, min_arc = _next_arc;

        if (_curr_length > 0 && _minor_count < _minor_limit) {

          // Minor iteration: select the best eligible arc from the

          // current candidate list

          ++_minor_count;

          min = 0;

              min = c;

              min_arc = e;

            }

            if (c >= 0) {

              _candidates[i--] = _candidates[--_curr_length];

            }

          }

          if (min < 0) {

            _in_arc = min_arc;

            return true;

          }

        }

        // Major iteration: build a new candidate list

        min = 0;

        _curr_length = 0;

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < 0) {

            _candidates[_curr_length++] = e;

            if (c < min) {

              min = c;

              min_arc = e;

            }

            if (_curr_length == _list_length) break;

          }

        }

        if (_curr_length < _list_length) {

          for (e = 0; e < _next_arc; ++e) {

            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

            if (c < 0) {

              _candidates[_curr_length++] = e;

              if (c < min) {

    }; //class CandidateListPivotRule

    // Implementation of the Altering Candidate List pivot rule

    class AlteringListPivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const IntVector  &_state;

      const CostVector &_pi;

      int &_in_arc;

      // Pivot rule data

      int _block_size, _head_length, _curr_length;

      int _next_arc;

      IntVector _candidates;

      CostVector _cand_cost;

      // Functor class to compare arcs during sort of the candidate list

      class SortFunc

      {

      private:

        const CostVector &_map;

      public:

        SortFunc(const CostVector &map) : _map(map) {}

        bool operator()(int left, int right) {

          return _map[left] > _map[right];

        }

      };

      SortFunc _sort_func;

    public:

      // Constructor

      AlteringListPivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

      {

        // The main parameters of the pivot rule

        const double BLOCK_SIZE_FACTOR = 1.5;

        const int MIN_BLOCK_SIZE = 10;

        const double HEAD_LENGTH_FACTOR = 0.1;

        const int MIN_HEAD_LENGTH = 3;

        _block_size = std::max( int(BLOCK_SIZE_FACTOR *

                                MIN_BLOCK_SIZE );

        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),

                                 MIN_HEAD_LENGTH );

        _candidates.resize(_head_length + _block_size);

        _curr_length = 0;

      }

      // Find next entering arc

      bool findEnteringArc() {

        // Check the current candidate list

        int e;

        for (int i = 0; i < _curr_length; ++i) {

          e = _candidates[i];

          _cand_cost[e] = _state[e] *

            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (_cand_cost[e] >= 0) {

            _candidates[i--] = _candidates[--_curr_length];

          }

        }

        // Extend the list

        int cnt = _block_size;

        int last_arc = 0;

        int limit = _head_length;

          _cand_cost[e] = _state[e] *

            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (_cand_cost[e] < 0) {

            _candidates[_curr_length++] = e;

            last_arc = e;

          }

          if (--cnt == 0) {

            if (_curr_length > limit) break;

            limit = 0;

            cnt = _block_size;

          }

        }

        if (_curr_length <= limit) {

          for (int e = 0; e < _next_arc; ++e) {

            _cand_cost[e] = _state[e] *

              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

    NetworkSimplex(const GR& graph) :

      _graph(graph), _node_id(graph), _arc_id(graph),

      INF(std::numeric_limits<Value>::has_infinity ?

          std::numeric_limits<Value>::infinity() :

          std::numeric_limits<Value>::max())

    {

      // Check the value types

      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,

        "The flow type of NetworkSimplex must be signed");

      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,

        "The cost type of NetworkSimplex must be signed");

      // Resize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      int all_node_num = _node_num + 1;

      _supply.resize(all_node_num);

      _pi.resize(all_node_num);

      _parent.resize(all_node_num);

      _pred.resize(all_node_num);

      _forward.resize(all_node_num);

      _thread.resize(all_node_num);

      _rev_thread.resize(all_node_num);

      _succ_num.resize(all_node_num);

      _last_succ.resize(all_node_num);

      // Copy the graph (store the arcs in a mixed order)

      int i = 0;

      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

        _node_id[n] = i;

      }

      int k = std::max(int(std::sqrt(double(_arc_num))), 10);

      i = 0;

      for (ArcIt a(_graph); a != INVALID; ++a) {

        _arc_id[a] = i;

        _source[i] = _node_id[_graph.source(a)];

        _target[i] = _node_id[_graph.target(a)];

        if ((i += k) >= _arc_num) i = (i % k) + 1;

      }

      // Initialize maps

            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;

          } else {

            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;

          }

          _supply[_source[i]] -= c;

          _supply[_target[i]] += c;

        }

      } else {

        for (int i = 0; i != _arc_num; ++i) {

          _cap[i] = _upper[i];

        }

      }

      // Initialize artifical cost

      Cost ART_COST;

      if (std::numeric_limits<Cost>::is_exact) {

      } else {

        ART_COST = std::numeric_limits<Cost>::min();

        for (int i = 0; i != _arc_num; ++i) {

          if (_cost[i] > ART_COST) ART_COST = _cost[i];

        }

        ART_COST = (ART_COST + 1) * _node_num;

      }

      // Initialize arc maps

      for (int i = 0; i != _arc_num; ++i) {

        _flow[i] = 0;

        _state[i] = STATE_LOWER;

      }

      // Set data for the artificial root node

      _root = _node_num;

      _parent[_root] = -1;

      _pred[_root] = -1;

      _thread[_root] = 0;

      _rev_thread[0] = _root;

      _succ_num[_root] = _node_num + 1;

      _last_succ[_root] = _root - 1;

      _supply[_root] = -_sum_supply;

      // Add artificial arcs and initialize the spanning tree data structure

      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

        _parent[u] = _root;

        _pred[u] = e;

        _thread[u] = u + 1;

        _rev_thread[u + 1] = u;

        _succ_num[u] = 1;

        _last_succ[u] = u;

        _cap[e] = INF;

        _state[e] = STATE_TREE;

          _flow[e] = _supply[u];

          _forward[u] = true;

        } else {

          _flow[e] = -_supply[u];

        }

      }

      return true;

    }

    // Find the join node

    void findJoinNode() {

      int u = _source[in_arc];

      int v = _target[in_arc];

      while (u != v) {

        if (_succ_num[u] < _succ_num[v]) {

          u = _parent[u];

        } else {

          v = _parent[v];

        }

      }

      join = u;

    ProblemType start() {

      PivotRuleImpl pivot(*this);