Improvments in min cost flow algorithms
authordeba
Fri, 15 Jun 2007 14:36:24 +0000
changeset 24578c791ee69a45
parent 2456 717a5134ddeb
child 2458 93b4132ac1e8
Improvments in min cost flow algorithms
- improved cycle cancelling
lemon/capacity_scaling.h
lemon/cycle_canceling.h
lemon/network_simplex.h
     1.1 --- a/lemon/capacity_scaling.h	Fri Jun 15 14:32:48 2007 +0000
     1.2 +++ b/lemon/capacity_scaling.h	Fri Jun 15 14:36:24 2007 +0000
     1.3 @@ -31,13 +31,15 @@
     1.4  
     1.5  #define WITH_SCALING
     1.6  
     1.7 +//#define _DEBUG_ITER_
     1.8 +
     1.9  namespace lemon {
    1.10  
    1.11    /// \addtogroup min_cost_flow
    1.12    /// @{
    1.13  
    1.14    /// \brief Implementation of the capacity scaling version of the
    1.15 -  /// succesive shortest path algorithm for finding a minimum cost flow.
    1.16 +  /// successive shortest path algorithm for finding a minimum cost flow.
    1.17    ///
    1.18    /// \ref lemon::CapacityScaling "CapacityScaling" implements a
    1.19    /// capacity scaling algorithm for finding a minimum cost flow.
    1.20 @@ -299,6 +301,9 @@
    1.21      /// \brief Map for identifing deficit nodes.
    1.22      DeficitBoolMap delta_deficit;
    1.23  
    1.24 +    /// \brief The deficit nodes.
    1.25 +    std::vector<ResNode> deficit_nodes;
    1.26 +
    1.27  #else //WITHOUT_CAPACITY_SCALING
    1.28      typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
    1.29        ResDijkstra;
    1.30 @@ -510,9 +515,9 @@
    1.31        return c;
    1.32      }
    1.33  
    1.34 -    /// \brief Runs the successive shortest path algorithm.
    1.35 +    /// \brief Runs the algorithm.
    1.36      ///
    1.37 -    /// Runs the successive shortest path algorithm.
    1.38 +    /// Runs the algorithm.
    1.39      ///
    1.40      /// \return \c true if a feasible flow can be found.
    1.41      bool run() {
    1.42 @@ -531,11 +536,13 @@
    1.43  
    1.44  #ifdef WITH_SCALING
    1.45        // Initilaizing delta value
    1.46 -      Capacity max_cap = 0;
    1.47 -      for (EdgeIt e(graph); e != INVALID; ++e) {
    1.48 -	if (capacity[e] > max_cap) max_cap = capacity[e];
    1.49 +      Supply max_sup = 0, max_dem = 0;
    1.50 +      for (NodeIt n(graph); n != INVALID; ++n) {
    1.51 +	if (supply[n] >  max_sup) max_sup =  supply[n];
    1.52 +	if (supply[n] < -max_dem) max_dem = -supply[n];
    1.53        }
    1.54 -      for (delta = 1; 2 * delta < max_cap; delta *= 2) ;
    1.55 +      if (max_dem < max_sup) max_sup = max_dem;
    1.56 +      for (delta = 1; 2 * delta < max_sup; delta *= 2) ;
    1.57  #endif
    1.58        return true;
    1.59      }
    1.60 @@ -546,6 +553,9 @@
    1.61      bool start() {
    1.62        typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
    1.63        typedef typename DeltaResGraph::Edge DeltaResEdge;
    1.64 +#ifdef _DEBUG_ITER_
    1.65 +      int dijk_num = 0;
    1.66 +#endif
    1.67  
    1.68        // Processing capacity scaling phases
    1.69        ResNode s, t;
    1.70 @@ -564,18 +574,35 @@
    1.71  
    1.72  	// Finding excess nodes
    1.73  	excess_nodes.clear();
    1.74 +	deficit_nodes.clear();
    1.75  	for (ResNodeIt n(res_graph); n != INVALID; ++n) {
    1.76 -	  if (imbalance[n] >= delta) excess_nodes.push_back(n);
    1.77 +	  if (imbalance[n] >=  delta) excess_nodes.push_back(n);
    1.78 +	  if (imbalance[n] <= -delta) deficit_nodes.push_back(n);
    1.79  	}
    1.80  	next_node = 0;
    1.81  
    1.82 -	// Finding successive shortest paths
    1.83 +	// Finding shortest paths
    1.84  	while (next_node < excess_nodes.size()) {
    1.85 +	  // Checking deficit nodes
    1.86 +	  if (delta > 1) {
    1.87 +	    bool delta_def = false;
    1.88 +	    for (int i = 0; i < deficit_nodes.size(); ++i) {
    1.89 +	      if (imbalance[deficit_nodes[i]] <= -delta) {
    1.90 +		delta_def = true;
    1.91 +		break;
    1.92 +	      }
    1.93 +	    }
    1.94 +	    if (!delta_def) break;
    1.95 +	  }
    1.96 +
    1.97  	  // Running Dijkstra
    1.98  	  s = excess_nodes[next_node];
    1.99  	  updater.init();
   1.100  	  dijkstra.init();
   1.101  	  dijkstra.addSource(s);
   1.102 +#ifdef _DEBUG_ITER_
   1.103 +	  ++dijk_num;
   1.104 +#endif
   1.105  	  if ((t = dijkstra.start(delta_deficit)) == INVALID) {
   1.106  	    if (delta > 1) {
   1.107  	      ++next_node;
   1.108 @@ -608,6 +635,10 @@
   1.109  	  if (imbalance[s] < delta) ++next_node;
   1.110  	}
   1.111        }
   1.112 +#ifdef _DEBUG_ITER_
   1.113 +      std::cout << "Cost Scaling algorithm finished with running Dijkstra algorithm "
   1.114 +	<< dijk_num << " times." << std::endl;
   1.115 +#endif
   1.116  
   1.117        // Handling nonzero lower bounds
   1.118        if (lower) {
   1.119 @@ -628,7 +659,7 @@
   1.120        if (excess_nodes.size() == 0) return true;
   1.121        next_node = 0;
   1.122  
   1.123 -      // Finding successive shortest paths
   1.124 +      // Finding shortest paths
   1.125        ResNode s, t;
   1.126        while ( imbalance[excess_nodes[next_node]] > 0 ||
   1.127  	      ++next_node < excess_nodes.size() )
     2.1 --- a/lemon/cycle_canceling.h	Fri Jun 15 14:32:48 2007 +0000
     2.2 +++ b/lemon/cycle_canceling.h	Fri Jun 15 14:36:24 2007 +0000
     2.3 @@ -22,13 +22,13 @@
     2.4  /// \ingroup min_cost_flow
     2.5  ///
     2.6  /// \file
     2.7 -/// \brief A cycle canceling algorithm for finding a minimum cost flow.
     2.8 +/// \brief A cycle-canceling algorithm for finding a minimum cost flow.
     2.9  
    2.10  #include <vector>
    2.11  #include <lemon/circulation.h>
    2.12  #include <lemon/graph_adaptor.h>
    2.13  
    2.14 -/// \brief The used cycle canceling method.
    2.15 +/// \brief The used cycle-canceling method.
    2.16  #define LIMITED_CYCLE_CANCELING
    2.17  //#define MIN_MEAN_CYCLE_CANCELING
    2.18  
    2.19 @@ -36,17 +36,21 @@
    2.20    #include <lemon/bellman_ford.h>
    2.21    /// \brief The maximum number of iterations for the first execution
    2.22    /// of the \ref lemon::BellmanFord "Bellman-Ford" algorithm.
    2.23 +  /// It should be at least 2.
    2.24    #define STARTING_LIMIT	2
    2.25    /// \brief The iteration limit for the
    2.26    /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
    2.27 -  /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
    2.28 +  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> in every round.
    2.29 +  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> must be greater than 1.
    2.30    #define ALPHA_MUL		3
    2.31    /// \brief The iteration limit for the
    2.32    /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
    2.33 -  /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
    2.34 +  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> in every round.
    2.35 +  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> must be greater than 1.
    2.36    #define ALPHA_DIV		2
    2.37  
    2.38  //#define _ONLY_ONE_CYCLE_
    2.39 +//#define _NO_BACK_STEP_
    2.40  //#define _DEBUG_ITER_
    2.41  #endif
    2.42  
    2.43 @@ -60,11 +64,11 @@
    2.44    /// \addtogroup min_cost_flow
    2.45    /// @{
    2.46  
    2.47 -  /// \brief Implementation of a cycle canceling algorithm for finding
    2.48 +  /// \brief Implementation of a cycle-canceling algorithm for finding
    2.49    /// a minimum cost flow.
    2.50    ///
    2.51 -  /// \ref lemon::CycleCanceling "CycleCanceling" implements a cycle
    2.52 -  /// canceling algorithm for finding a minimum cost flow.
    2.53 +  /// \ref lemon::CycleCanceling "CycleCanceling" implements a
    2.54 +  /// cycle-canceling algorithm for finding a minimum cost flow.
    2.55    ///
    2.56    /// \param Graph The directed graph type the algorithm runs on.
    2.57    /// \param LowerMap The type of the lower bound map.
    2.58 @@ -118,26 +122,6 @@
    2.59  
    2.60    protected:
    2.61  
    2.62 -    /// \brief Map adaptor class for demand map.
    2.63 -    class DemandMap : public MapBase<Node, Supply>
    2.64 -    {
    2.65 -    private:
    2.66 -
    2.67 -      const SupplyMap *map;
    2.68 -
    2.69 -    public:
    2.70 -
    2.71 -      typedef typename MapBase<Node, Supply>::Value Value;
    2.72 -      typedef typename MapBase<Node, Supply>::Key Key;
    2.73 -
    2.74 -      DemandMap(const SupplyMap &_map) : map(&_map) {}
    2.75 -
    2.76 -      Value operator[](const Key &e) const {
    2.77 -	return -(*map)[e];
    2.78 -      }
    2.79 -
    2.80 -    }; //class DemandMap
    2.81 -
    2.82      /// \brief Map adaptor class for handling residual edge costs.
    2.83      class ResCostMap : public MapBase<ResEdge, Cost>
    2.84      {
    2.85 @@ -170,8 +154,6 @@
    2.86      const CostMap &cost;
    2.87      /// \brief The modified supply map.
    2.88      SupplyRefMap supply;
    2.89 -    /// \brief The modified demand map.
    2.90 -    DemandMap demand;
    2.91      /// \brief The sum of supply values equals zero.
    2.92      bool valid_supply;
    2.93  
    2.94 @@ -199,7 +181,7 @@
    2.95  		    const CostMap &_cost,
    2.96  		    const SupplyMap &_supply ) :
    2.97        graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
    2.98 -      supply(_graph), demand(supply), flow(_graph, 0),
    2.99 +      supply(_graph), flow(_graph, 0),
   2.100        res_graph(_graph, capacity, flow), res_cost(cost)
   2.101      {
   2.102        // Removing nonzero lower bounds
   2.103 @@ -229,7 +211,7 @@
   2.104  		    const CostMap &_cost,
   2.105  		    const SupplyMap &_supply ) :
   2.106        graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
   2.107 -      supply(_supply), demand(supply), flow(_graph, 0),
   2.108 +      supply(_supply), flow(_graph, 0),
   2.109        res_graph(_graph, capacity, flow), res_cost(cost)
   2.110      {
   2.111        // Checking the sum of supply values
   2.112 @@ -259,7 +241,7 @@
   2.113  		    Node _s, Node _t,
   2.114  		    Supply _flow_value ) :
   2.115        graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   2.116 -      supply(_graph), demand(supply), flow(_graph, 0),
   2.117 +      supply(_graph), flow(_graph, 0),
   2.118        res_graph(_graph, capacity, flow), res_cost(cost)
   2.119      {
   2.120        // Removing nonzero lower bounds
   2.121 @@ -295,7 +277,7 @@
   2.122  		    Node _s, Node _t,
   2.123  		    Supply _flow_value ) :
   2.124        graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
   2.125 -      supply(_graph, 0), demand(supply), flow(_graph, 0),
   2.126 +      supply(_graph, 0), flow(_graph, 0),
   2.127        res_graph(_graph, capacity, flow), res_cost(cost)
   2.128      {
   2.129        supply[_s] =  _flow_value;
   2.130 @@ -345,14 +327,14 @@
   2.131  
   2.132        // Finding a feasible flow
   2.133        Circulation< Graph, Capacity, FlowMap, ConstMap<Edge, Capacity>,
   2.134 -		   CapacityRefMap, DemandMap >
   2.135 +		   CapacityRefMap, SupplyMap >
   2.136  	circulation( graph, constMap<Edge>((Capacity)0),
   2.137 -		     capacity, demand, flow );
   2.138 +		     capacity, supply, flow );
   2.139        return circulation.run() == -1;
   2.140      }
   2.141  
   2.142  #ifdef LIMITED_CYCLE_CANCELING
   2.143 -    /// \brief Executes a cycle canceling algorithm using
   2.144 +    /// \brief Executes a cycle-canceling algorithm using
   2.145      /// \ref lemon::BellmanFord "Bellman-Ford" algorithm with limited
   2.146      /// iteration count.
   2.147      bool start() {
   2.148 @@ -373,8 +355,13 @@
   2.149  	int iter_num = 0;
   2.150  	bool cycle_found = false;
   2.151  	while (!cycle_found) {
   2.152 +#ifdef _NO_BACK_STEP_
   2.153 +	  int curr_iter_num = length_bound <= node_num ?
   2.154 +			      length_bound - iter_num : node_num - iter_num;
   2.155 +#else
   2.156  	  int curr_iter_num = iter_num + length_bound <= node_num ?
   2.157  			      length_bound : node_num - iter_num;
   2.158 +#endif
   2.159  	  iter_num += curr_iter_num;
   2.160  	  int real_iter_num = curr_iter_num;
   2.161  	  for (int i = 0; i < curr_iter_num; ++i) {
   2.162 @@ -432,7 +419,7 @@
   2.163        }
   2.164  
   2.165  #ifdef _DEBUG_ITER_
   2.166 -      std::cout << "Limited cycle canceling algorithm finished. "
   2.167 +      std::cout << "Limited cycle-canceling algorithm finished. "
   2.168  		<< "Found " << cycle_num << " negative cycles."
   2.169  		<< std::endl;
   2.170  #endif
   2.171 @@ -447,7 +434,7 @@
   2.172  #endif
   2.173  
   2.174  #ifdef MIN_MEAN_CYCLE_CANCELING
   2.175 -    /// \brief Executes the minimum mean cycle canceling algorithm
   2.176 +    /// \brief Executes the minimum mean cycle-canceling algorithm
   2.177      /// using \ref lemon::MinMeanCycle "MinMeanCycle" class.
   2.178      bool start() {
   2.179        typedef Path<ResGraph> ResPath;
   2.180 @@ -486,7 +473,7 @@
   2.181        }
   2.182  
   2.183  #ifdef _DEBUG_ITER_
   2.184 -      std::cout << "Minimum mean cycle canceling algorithm finished. "
   2.185 +      std::cout << "Minimum mean cycle-canceling algorithm finished. "
   2.186  		<< "Found " << cycle_num << " negative cycles."
   2.187  		<< std::endl;
   2.188  #endif
     3.1 --- a/lemon/network_simplex.h	Fri Jun 15 14:32:48 2007 +0000
     3.2 +++ b/lemon/network_simplex.h	Fri Jun 15 14:36:24 2007 +0000
     3.3 @@ -275,6 +275,8 @@
     3.4        if (!(valid_supply = sum == 0)) return;
     3.5  
     3.6        // Copying graph_ref to graph
     3.7 +      graph.reserveNode(countNodes(graph_ref) + 1);
     3.8 +      graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));
     3.9        copyGraph(graph, graph_ref)
    3.10  	.edgeMap(cost, _cost)
    3.11  	.nodeRef(node_ref)
    3.12 @@ -477,7 +479,9 @@
    3.13        root = graph.addNode();
    3.14        parent[root] = INVALID;
    3.15        pred_edge[root] = INVALID;
    3.16 -      depth[root] = supply[root] = potential[root] = 0;
    3.17 +      depth[root] = 0;
    3.18 +      supply[root] = 0;
    3.19 +      potential[root] = 0;
    3.20  
    3.21        // Adding artificial edges to the graph and initializing the node
    3.22        // maps of the spanning tree data structure
    3.23 @@ -513,7 +517,7 @@
    3.24  #ifdef EDGE_BLOCK_PIVOT
    3.25        // Initializing block_size for the edge block pivot rule
    3.26        int edge_num = countEdges(graph);
    3.27 -      block_size = edge_num >= BLOCK_NUM * MIN_BLOCK_SIZE ? 
    3.28 +      block_size = edge_num >= BLOCK_NUM * MIN_BLOCK_SIZE ?
    3.29  		   edge_num / BLOCK_NUM : MIN_BLOCK_SIZE;
    3.30  #endif
    3.31  #ifdef CANDIDATE_LIST_PIVOT