RhodeCode

@defgroup algs Algorithms

\brief This group describes the several algorithms

implemented in LEMON.

This group describes the several algorithms

implemented in LEMON.

*/

/**

@defgroup search Graph Search

@ingroup algs

\brief Common graph search algorithms.

This group describes the common graph search algorithms, namely

\e breadth-first \e search (BFS) and \e depth-first \e search (DFS).

*/

/**

@defgroup shortest_path Shortest Path Algorithms

@ingroup algs

\brief Algorithms for finding shortest paths.

This group describes the algorithms for finding shortest paths in digraphs.

 - \ref Dijkstra algorithm for finding shortest paths from a source node

   when all arc lengths are non-negative.

 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths

   from a source node when arc lenghts can be either positive or negative,

   but the digraph should not contain directed cycles with negative total

   length.

 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms

   for solving the \e all-pairs \e shortest \e paths \e problem when arc

   lenghts can be either positive or negative, but the digraph should

   not contain directed cycles with negative total length.

 - \ref Suurballe A successive shortest path algorithm for finding

   arc-disjoint paths between two nodes having minimum total length.

*/

/**

@defgroup max_flow Maximum Flow Algorithms

@ingroup algs

\brief Algorithms for finding maximum flows.

This group describes the algorithms for finding maximum flows and

feasible circulations.

The \e maximum \e flow \e problem is to find a flow of maximum value between

a single source and a single target. Formally, there is a \f$G=(V,A)\f$

\f$s, t \in V\f$ source and target nodes.

following optimization problem.

LEMON contains several algorithms for solving maximum flow problems:

- \ref EdmondsKarp Edmonds-Karp algorithm.

- \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

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

problem of the maximum flow.

*/

/**

@defgroup min_cost_flow Minimum Cost Flow Algorithms

@ingroup algs

\brief Algorithms for finding minimum cost flows and circulations.

circulations.

The \e minimum \e cost \e flow \e problem is to find a feasible flow of

minimum total cost from a set of supply nodes to a set of demand nodes

Formally, let \f$G=(V,A)\f$ be a digraph,

\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and

\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow

   capacity scaling.

*/

/**

@defgroup min_cut Minimum Cut Algorithms

@ingroup algs

\brief Algorithms for finding minimum cut in graphs.

This group describes 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

outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a

\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum

cut is the \f$X\f$ solution of the next optimization problem:

\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}

    \sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f]

LEMON contains several algorithms related to minimum cut problems:

- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut

  in directed graphs.

- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for

  calculating minimum cut in undirected graphs.

- \ref GomoryHuTree "Gomory-Hu tree computation" for calculating

  all-pairs minimum cut in undirected graphs.

If you want to find minimum cut just between two distinict nodes,

see the \ref max_flow "maximum flow problem".

*/

/**

@defgroup graph_prop Connectivity and Other Graph Properties

@ingroup algs

\brief Algorithms for discovering the graph properties

This group describes the algorithms for discovering the graph properties

like connectivity, bipartiteness, euler property, simplicity etc.

\image html edge_biconnected_components.png

\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth

*/

/**

@defgroup planar Planarity Embedding and Drawing

@ingroup algs

\brief Algorithms for planarity checking, embedding and drawing

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Copyright (C) 2003-2009

 * 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

/// \ingroup min_cost_flow

///

/// \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 {

  /// \addtogroup min_cost_flow

  /// @{

  /// \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.

  ///

  /// \tparam GR The digraph type the algorithm runs on.

  /// \tparam F The value type used for flow amounts, capacity bounds

  /// and supply values in the algorithm. By default it is \c int.

  /// \tparam C The value type used for costs and potentials in the

  /// algorithm. By default it is the same as \c F.

  ///

  /// \warning Both value types must be signed and all input data must

  /// be integer.

  ///

  /// \note %NetworkSimplex provides five different pivot rule

  template <typename GR, typename F = int, typename C = F>

  class NetworkSimplex

  {

  public:

    /// The flow type of the algorithm

    typedef F Flow;

    /// The cost type of the algorithm

    typedef C Cost;

    /// The type of the flow map

    typedef typename GR::template ArcMap<Flow> FlowMap;

    /// The type of the potential map

    typedef typename GR::template NodeMap<Cost> PotentialMap;

  public:

    /// \brief Enum type for selecting the pivot rule.

    ///

    /// Enum type 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 {

      /// The First Eligible pivot rule.

      /// The next eligible arc is selected in a wraparound fashion

      /// in every iteration.

      FIRST_ELIGIBLE,

      /// The Best Eligible pivot rule.

      /// The best eligible arc is selected in every iteration.

      BEST_ELIGIBLE,

      /// The Block Search pivot rule.

      /// A specified number of arcs are examined in every iteration

      /// in a wraparound fashion and the best eligible arc is selected

      /// from this block.

      BLOCK_SEARCH,

      /// The Candidate List pivot rule.

      /// In a major iteration a candidate list is built from eligible arcs

      /// in a wraparound fashion and in the following minor iterations

      /// the best eligible arc is selected from this list.

      CANDIDATE_LIST,

      /// The Altering Candidate List pivot rule.

      /// It is a modified version of the Candidate List method.

      /// It keeps only the several best eligible arcs from the former

      /// candidate list and extends this list in every iteration.

      ALTERING_LIST

    };

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(GR);

    typedef typename GR::template ArcMap<Flow> FlowArcMap;

    typedef typename GR::template ArcMap<Cost> CostArcMap;

    typedef typename GR::template NodeMap<Flow> FlowNodeMap;

    typedef std::vector<Arc> ArcVector;

    typedef std::vector<Node> NodeVector;

    typedef std::vector<int> IntVector;

    typedef std::vector<bool> BoolVector;

    typedef std::vector<Flow> FlowVector;

    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

    FlowArcMap *_plower;

    FlowArcMap *_pupper;

    CostArcMap *_pcost;

    FlowNodeMap *_psupply;

    bool _pstsup;

    Node _psource, _ptarget;

    Flow _pstflow;

    // Result maps

    FlowMap *_flow_map;

    PotentialMap *_potential_map;

    bool _local_flow;

    bool _local_potential;

    // Data structures for storing the digraph

    IntNodeMap _node_id;

    ArcVector _arc_ref;

    IntVector _source;

    IntVector _target;

    // Node and arc data

    FlowVector _cap;

    CostVector _cost;

    FlowVector _supply;

    FlowVector _flow;

    CostVector _pi;

    // Data for storing the spanning tree structure

    IntVector _parent;

    IntVector _pred;

    IntVector _thread;

    IntVector _rev_thread;

    IntVector _succ_num;

    IntVector _last_succ;

    IntVector _dirty_revs;

    BoolVector _forward;

    IntVector _state;

    int _root;

    // Temporary data used in the current pivot iteration

    int in_arc, join, u_in, v_in, u_out, v_out;

    int first, second, right, last;

    int stem, par_stem, new_stem;

    Flow delta;

  private:

    // Implementation of the First Eligible pivot rule

    class FirstEligiblePivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

          _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]]);

            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 == 0) return false;

        _next_arc = last_arc + 1;

        // Make heap of the candidate list (approximating a partial sort)

        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,

                   _sort_func );

        // Pop the first element of the heap

        _in_arc = _candidates[0];

        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,

                  _sort_func );

        _curr_length = std::min(_head_length, _curr_length - 1);

        return true;

      }

    }; //class AlteringListPivotRule

  public:

    /// \brief Constructor.

    ///

    ///

    /// \param graph The digraph the algorithm runs on.

    NetworkSimplex(const GR& graph) :

      _graph(graph),

      _plower(NULL), _pupper(NULL), _pcost(NULL),

      _flow_map(NULL), _potential_map(NULL),

      _local_flow(false), _local_potential(false),

      _node_id(graph)

    {

      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&

                   std::numeric_limits<Flow>::is_signed,

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

      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&

                   std::numeric_limits<Cost>::is_signed,

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

    }

    /// Destructor.

    ~NetworkSimplex() {

      if (_local_flow) delete _flow_map;

      if (_local_potential) delete _potential_map;

    }

    /// \brief Set the lower bounds on the arcs.

    ///

    /// This function sets the lower bounds on the arcs.

    /// If neither this function nor \ref boundMaps() is used before

    /// calling \ref run(), the lower bounds will be set to zero

    /// on all arcs.

    ///

    /// \param map An arc map storing the lower bounds.

    /// Its \c Value type must be convertible to the \c Flow type

    /// of the algorithm.

    ///

    /// \return <tt>(*this)</tt>

    template <typename LOWER>

    NetworkSimplex& lowerMap(const LOWER& map) {

      delete _plower;

      _plower = new FlowArcMap(_graph);

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

        (*_plower)[a] = map[a];

      }

      return *this;

    }

    /// \brief Set the upper bounds (capacities) on the arcs.

    ///

    /// This function sets the upper bounds (capacities) on the arcs.

    /// If none of the functions \ref upperMap(), \ref capacityMap()

    /// and \ref boundMaps() is used before calling \ref run(),

    /// the upper bounds (capacities) will be set to

    /// \c std::numeric_limits<Flow>::max() on all arcs.

    ///

    /// \param map An arc map storing the upper bounds.

    /// Its \c Value type must be convertible to the \c Flow type

    /// of the algorithm.

    ///

    /// \return <tt>(*this)</tt>

    template<typename UPPER>

    NetworkSimplex& upperMap(const UPPER& map) {

      delete _pupper;

      _pupper = new FlowArcMap(_graph);

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

        (*_pupper)[a] = map[a];

      }

      return *this;

    }

    /// \brief Set the upper bounds (capacities) on the arcs.

    ///

    /// This function sets the upper bounds (capacities) on the arcs.

    }

    /// \brief Set the supply values of the nodes.

    ///

    /// This function sets the supply values of the nodes.

    /// If neither this function nor \ref stSupply() is used before

    /// calling \ref run(), the supply of each node will be set to zero.

    /// (It makes sense only if non-zero lower bounds are given.)

    ///

    /// \param map A node map storing the supply values.

    /// Its \c Value type must be convertible to the \c Flow type

    /// of the algorithm.

    ///

    /// \return <tt>(*this)</tt>

    template<typename SUP>

    NetworkSimplex& supplyMap(const SUP& map) {

      delete _psupply;

      _pstsup = false;

      _psupply = new FlowNodeMap(_graph);

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

        (*_psupply)[n] = map[n];

      }

      return *this;

    }

    /// \brief Set single source and target nodes and a supply value.

    ///

    /// This function sets a single source node and a single target node

    /// and the required flow value.

    /// If neither this function nor \ref supplyMap() is used before

    /// calling \ref run(), the supply of each node will be set to zero.

    /// (It makes sense only if non-zero lower bounds are given.)

    ///

    /// \param s The source node.

    /// \param t The target node.

    /// \param k The required amount of flow from node \c s to node \c t

    /// (i.e. the supply of \c s and the demand of \c t).

    ///

    /// \return <tt>(*this)</tt>

    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {

      delete _psupply;

      _psupply = NULL;

      _pstsup = true;

      _psource = s;

      _ptarget = t;

      _pstflow = k;

      return *this;

    }

    /// \brief Set the flow map.

    ///

    /// This function sets the flow map.

    /// If it is not used before calling \ref run(), an instance will

    /// be allocated automatically. The destructor deallocates this

    /// automatically allocated map, of course.

    ///

    /// \return <tt>(*this)</tt>

    NetworkSimplex& flowMap(FlowMap& map) {

      if (_local_flow) {

        delete _flow_map;

        _local_flow = false;

      }

      _flow_map = ↦

      return *this;

    }

    /// \brief Set the potential map.

    ///

    /// This function sets the potential map, which is used for storing

    /// the dual solution.

    /// If it is not used before calling \ref run(), an instance will

    /// be allocated automatically. The destructor deallocates this

    /// automatically allocated map, of course.

    ///

    /// \return <tt>(*this)</tt>

    NetworkSimplex& potentialMap(PotentialMap& map) {

      if (_local_potential) {

        delete _potential_map;

        _local_potential = false;

      }

      _potential_map = ↦

      return *this;

    }

    /// \name Execution Control

    /// The algorithm can be executed using \ref run().

    /// @{

    /// \brief Run the algorithm.

    ///

    /// This function runs the algorithm.

    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),

    /// \code

    ///   NetworkSimplex<ListDigraph> ns(graph);

    ///   ns.boundMaps(lower, upper).costMap(cost)

    ///     .supplyMap(sup).run();

    /// \endcode

    ///

    /// This function can be called more than once. All the parameters

    /// that have been given are kept for the next call, unless

    /// \ref reset() is called, thus only the modified parameters

    /// have to be set again. See \ref reset() for examples.

    ///

    /// \param pivot_rule The pivot rule that will be used during the

    /// algorithm. For more information see \ref PivotRule.

    ///

    /// \return \c true if a feasible flow can be found.

    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {

      return init() && start(pivot_rule);

    }

    /// \brief Reset all the parameters that have been given before.

    ///

    /// This function resets all the paramaters that have been given

    ///

    /// It is useful for multiple run() calls. If this function is not

    /// used, all the parameters given before are kept for the next

    /// \ref run() call.

    ///

    /// For example,

    /// \code

    ///   NetworkSimplex<ListDigraph> ns(graph);

    ///

    ///   // First run

    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)

    ///     .supplyMap(sup).run();

    ///

    ///   // Run again with modified cost map (reset() is not called,

    ///   // so only the cost map have to be set again)

    ///   cost[e] += 100;

    ///   ns.costMap(cost).run();

    ///

    ///   // Run again from scratch using reset()

    ///   // (the lower bounds will be set to zero on all arcs)

    ///   ns.reset();

    ///   ns.capacityMap(cap).costMap(cost)

    ///     .supplyMap(sup).run();

    /// \endcode

    ///

    /// \return <tt>(*this)</tt>

    NetworkSimplex& reset() {

      delete _plower;

      delete _pupper;

      delete _pcost;

      delete _psupply;

      _plower = NULL;

      _pupper = NULL;

      _pcost = NULL;

      _psupply = NULL;

      _pstsup = false;

      return *this;

    }

    /// @}

    /// \name Query Functions

    /// The results of the algorithm can be obtained using these

    /// functions.\n

    /// The \ref run() function must be called before using them.

    /// @{

    /// \brief Return the total cost of the found flow.

    ///

    /// This function returns the total cost of the found flow.

    /// The complexity of the function is O(e).

    ///

    /// \note The return type of the function can be specified as a

    /// template parameter. For example,

    /// \code

    ///   ns.totalCost<double>();

    /// \endcode

    /// It is useful if the total cost cannot be stored in the \c Cost

    /// type of the algorithm, which is the default return type of the

    /// function.

    ///

    /// \pre \ref run() must be called before using this function.

    template <typename Num>

    Num totalCost() const {

      Num c = 0;

      if (_pcost) {

        for (ArcIt e(_graph); e != INVALID; ++e)

          c += (*_flow_map)[e] * (*_pcost)[e];

      } else {

        for (ArcIt e(_graph); e != INVALID; ++e)

          c += (*_flow_map)[e];

      }

      return c;

    }

#ifndef DOXYGEN

    Cost totalCost() const {

      return totalCost<Cost>();

    }

#endif

    /// \brief Return the flow on the given arc.

    ///

    const PotentialMap& potentialMap() const {

      return *_potential_map;

    }

    /// @}

  private:

    // Initialize internal data structures

    bool init() {

      // Initialize result maps

      if (!_flow_map) {

        _flow_map = new FlowMap(_graph);

        _local_flow = true;

      }

      if (!_potential_map) {

        _potential_map = new PotentialMap(_graph);

        _local_potential = true;

      }

      // Initialize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      int all_node_num = _node_num + 1;

      int all_arc_num = _arc_num + _node_num;

      if (_node_num == 0) return false;

      _arc_ref.resize(_arc_num);

      _source.resize(all_arc_num);

      _target.resize(all_arc_num);

      _cap.resize(all_arc_num);

      _cost.resize(all_arc_num);

      _supply.resize(all_node_num);

      _flow.resize(all_arc_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);

      _state.resize(all_arc_num);

      // Initialize node related data

      bool valid_supply = true;

      if (!_pstsup && !_psupply) {

        _pstsup = true;

        _psource = _ptarget = NodeIt(_graph);

        _pstflow = 0;

      }

      if (_psupply) {

        int i = 0;

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

          _node_id[n] = i;

          _supply[i] = (*_psupply)[n];

        }

      } else {

        int i = 0;

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

          _node_id[n] = i;

          _supply[i] = 0;

        }

        _supply[_node_id[_psource]] =  _pstflow;

      }

      if (!valid_supply) return false;

      // 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] = all_node_num;

      _last_succ[_root] = _root - 1;

      // Store the arcs in a mixed order

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

      int i = 0;

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

        _arc_ref[i] = e;

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

      }

      // Initialize arc maps

      if (_pupper && _pcost) {

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

          Arc e = _arc_ref[i];

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

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

          _cap[i] = (*_pupper)[e];

          _cost[i] = (*_pcost)[e];

          _flow[i] = 0;

          _state[i] = STATE_LOWER;

        }

      } else {

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

          Arc e = _arc_ref[i];

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

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

          _flow[i] = 0;

          _state[i] = STATE_LOWER;

        }

        if (_pupper) {

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

            _cap[i] = (*_pupper)[_arc_ref[i]];

        } else {