lemon/capacity_scaling.h
changeset 2549 88b81ec599ed
parent 2533 aea952a1af99
child 2553 bfced05fa852
equal deleted inserted replaced
5:d3a88fa60e48 6:87226c005d1e
    23 ///
    23 ///
    24 /// \file
    24 /// \file
    25 /// \brief The capacity scaling algorithm for finding a minimum cost
    25 /// \brief The capacity scaling algorithm for finding a minimum cost
    26 /// flow.
    26 /// flow.
    27 
    27 
       
    28 #include <lemon/graph_adaptor.h>
       
    29 #include <lemon/bin_heap.h>
    28 #include <vector>
    30 #include <vector>
    29 #include <lemon/graph_adaptor.h>
       
    30 #include <lemon/dijkstra.h>
       
    31 
       
    32 #define WITH_SCALING
       
    33 
       
    34 #ifdef WITH_SCALING
       
    35 #define SCALING_FACTOR	  2
       
    36 #endif
       
    37 
       
    38 //#define _DEBUG_ITER_
       
    39 
    31 
    40 namespace lemon {
    32 namespace lemon {
    41 
    33 
    42   /// \addtogroup min_cost_flow
    34   /// \addtogroup min_cost_flow
    43   /// @{
    35   /// @{
    44 
    36 
    45   /// \brief Implementation of the capacity scaling version of the
    37   /// \brief Implementation of the capacity scaling version of the
    46   /// successive shortest path algorithm for finding a minimum cost
    38   /// successive shortest path algorithm for finding a minimum cost
    47   /// flow.
    39   /// flow.
    48   ///
    40   ///
    49   /// \ref lemon::CapacityScaling "CapacityScaling" implements the
    41   /// \ref CapacityScaling implements the capacity scaling version
    50   /// capacity scaling version of the successive shortest path
    42   /// of the successive shortest path algorithm for finding a minimum
    51   /// algorithm for finding a minimum cost flow.
    43   /// cost flow.
    52   ///
    44   ///
    53   /// \param Graph The directed graph type the algorithm runs on.
    45   /// \param Graph The directed graph type the algorithm runs on.
    54   /// \param LowerMap The type of the lower bound map.
    46   /// \param LowerMap The type of the lower bound map.
    55   /// \param CapacityMap The type of the capacity (upper bound) map.
    47   /// \param CapacityMap The type of the capacity (upper bound) map.
    56   /// \param CostMap The type of the cost (length) map.
    48   /// \param CostMap The type of the cost (length) map.
    57   /// \param SupplyMap The type of the supply map.
    49   /// \param SupplyMap The type of the supply map.
    58   ///
    50   ///
    59   /// \warning
    51   /// \warning
    60   /// - Edge capacities and costs should be nonnegative integers.
    52   /// - Edge capacities and costs should be nonnegative integers.
    61   ///	However \c CostMap::Value should be signed type.
    53   ///   However \c CostMap::Value should be signed type.
    62   /// - Supply values should be signed integers.
    54   /// - Supply values should be signed integers.
    63   /// - \c LowerMap::Value must be convertible to
    55   /// - \c LowerMap::Value must be convertible to
    64   ///	\c CapacityMap::Value and \c CapacityMap::Value must be
    56   ///   \c CapacityMap::Value and \c CapacityMap::Value must be
    65   ///	convertible to \c SupplyMap::Value.
    57   ///   convertible to \c SupplyMap::Value.
    66   ///
    58   ///
    67   /// \author Peter Kovacs
    59   /// \author Peter Kovacs
    68 
    60 
    69   template < typename Graph,
    61   template < typename Graph,
    70 	     typename LowerMap = typename Graph::template EdgeMap<int>,
    62              typename LowerMap = typename Graph::template EdgeMap<int>,
    71 	     typename CapacityMap = LowerMap,
    63              typename CapacityMap = LowerMap,
    72 	     typename CostMap = typename Graph::template EdgeMap<int>,
    64              typename CostMap = typename Graph::template EdgeMap<int>,
    73 	     typename SupplyMap = typename Graph::template NodeMap
    65              typename SupplyMap = typename Graph::template NodeMap
    74 				  <typename CapacityMap::Value> >
    66                                   <typename CapacityMap::Value> >
    75   class CapacityScaling
    67   class CapacityScaling
    76   {
    68   {
    77     typedef typename Graph::Node Node;
    69     typedef typename Graph::Node Node;
    78     typedef typename Graph::NodeIt NodeIt;
    70     typedef typename Graph::NodeIt NodeIt;
    79     typedef typename Graph::Edge Edge;
    71     typedef typename Graph::Edge Edge;
    85     typedef typename CapacityMap::Value Capacity;
    77     typedef typename CapacityMap::Value Capacity;
    86     typedef typename CostMap::Value Cost;
    78     typedef typename CostMap::Value Cost;
    87     typedef typename SupplyMap::Value Supply;
    79     typedef typename SupplyMap::Value Supply;
    88     typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
    80     typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
    89     typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
    81     typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
    90 
    82     typedef typename Graph::template NodeMap<Edge> PredMap;
    91     typedef ResGraphAdaptor< const Graph, Capacity,
       
    92 			     CapacityRefMap, CapacityRefMap > ResGraph;
       
    93     typedef typename ResGraph::Node ResNode;
       
    94     typedef typename ResGraph::NodeIt ResNodeIt;
       
    95     typedef typename ResGraph::Edge ResEdge;
       
    96     typedef typename ResGraph::EdgeIt ResEdgeIt;
       
    97 
    83 
    98   public:
    84   public:
       
    85 
       
    86     /// \brief Type to enable or disable capacity scaling.
       
    87     enum ScalingEnum {
       
    88       WITH_SCALING = 0,
       
    89       WITHOUT_SCALING = -1
       
    90     };
    99 
    91 
   100     /// \brief The type of the flow map.
    92     /// \brief The type of the flow map.
   101     typedef CapacityRefMap FlowMap;
    93     typedef CapacityRefMap FlowMap;
   102     /// \brief The type of the potential map.
    94     /// \brief The type of the potential map.
   103     typedef typename Graph::template NodeMap<Cost> PotentialMap;
    95     typedef typename Graph::template NodeMap<Cost> PotentialMap;
   104 
    96 
   105   protected:
    97   protected:
   106 
    98 
   107     /// \brief Map adaptor class for handling reduced edge costs.
    99     /// \brief Special implementation of the \ref Dijkstra algorithm
   108     class ReducedCostMap : public MapBase<ResEdge, Cost>
   100     /// for finding shortest paths in the residual network of the graph
       
   101     /// with respect to the reduced edge costs and modifying the
       
   102     /// node potentials according to the distance of the nodes.
       
   103     class ResidualDijkstra
   109     {
   104     {
   110     private:
   105       typedef typename Graph::template NodeMap<Cost> CostNodeMap;
   111 
   106       typedef typename Graph::template NodeMap<Edge> PredMap;
   112       const ResGraph &gr;
   107 
   113       const CostMap &cost_map;
   108       typedef typename Graph::template NodeMap<int> HeapCrossRef;
   114       const PotentialMap &pot_map;
   109       typedef BinHeap<Cost, HeapCrossRef> Heap;
       
   110 
       
   111     protected:
       
   112 
       
   113       /// \brief The directed graph the algorithm runs on.
       
   114       const Graph &graph;
       
   115 
       
   116       /// \brief The flow map.
       
   117       const FlowMap &flow;
       
   118       /// \brief The residual capacity map.
       
   119       const CapacityRefMap &res_cap;
       
   120       /// \brief The cost map.
       
   121       const CostMap &cost;
       
   122       /// \brief The excess map.
       
   123       const SupplyRefMap &excess;
       
   124 
       
   125       /// \brief The potential map.
       
   126       PotentialMap &potential;
       
   127 
       
   128       /// \brief The distance map.
       
   129       CostNodeMap dist;
       
   130       /// \brief The map of predecessors edges.
       
   131       PredMap &pred;
       
   132       /// \brief The processed (i.e. permanently labeled) nodes.
       
   133       std::vector<Node> proc_nodes;
   115 
   134 
   116     public:
   135     public:
   117 
   136 
   118       ReducedCostMap( const ResGraph &_gr,
   137       /// \brief The constructor of the class.
   119 		      const CostMap &_cost,
   138       ResidualDijkstra( const Graph &_graph,
   120 		      const PotentialMap &_pot ) :
   139                         const FlowMap &_flow,
   121 	gr(_gr), cost_map(_cost), pot_map(_pot) {}
   140                         const CapacityRefMap &_res_cap,
   122 
   141                         const CostMap &_cost,
   123       Cost operator[](const ResEdge &e) const {
   142                         const SupplyMap &_excess,
   124 	return ResGraph::forward(e) ?
   143                         PotentialMap &_potential,
   125 	   cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :
   144                         PredMap &_pred ) :
   126 	  -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
   145         graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost),
   127       }
   146         excess(_excess), potential(_potential), dist(_graph),
   128 
   147         pred(_pred)
   129     }; //class ReducedCostMap
   148       {}
   130 
   149 
   131     /// \brief Map class for the \ref lemon::Dijkstra "Dijkstra"
   150       /// \brief Runs the algorithm from the given source node.
   132     /// algorithm to update node potentials.
   151       Node run(Node s, Capacity delta) {
   133     class PotentialUpdateMap : public MapBase<ResNode, Cost>
   152         HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP);
   134     {
   153         Heap heap(heap_cross_ref);
   135     private:
   154         heap.push(s, 0);
   136 
   155         pred[s] = INVALID;
   137       PotentialMap *pot;
   156         proc_nodes.clear();
   138       typedef std::pair<ResNode, Cost> Pair;
   157 
   139       std::vector<Pair> data;
   158         // Processing nodes
   140 
   159         while (!heap.empty() && excess[heap.top()] > -delta) {
   141     public:
   160           Node u = heap.top(), v;
   142 
   161           Cost d = heap.prio() - potential[u], nd;
   143       void potentialMap(PotentialMap &_pot) {
   162           dist[u] = heap.prio();
   144 	pot = &_pot;
   163           heap.pop();
   145       }
   164           proc_nodes.push_back(u);
   146 
   165 
   147       void init() {
   166           // Traversing outgoing edges
   148 	data.clear();
   167           for (OutEdgeIt e(graph, u); e != INVALID; ++e) {
   149       }
   168             if (res_cap[e] >= delta) {
   150 
   169               v = graph.target(e);
   151       void set(const ResNode &n, const Cost &v) {
   170               switch(heap.state(v)) {
   152 	data.push_back(Pair(n, v));
   171               case Heap::PRE_HEAP:
   153       }
   172                 heap.push(v, d + cost[e] + potential[v]);
   154 
   173                 pred[v] = e;
   155       void update() {
   174                 break;
   156 	Cost val = data[data.size()-1].second;
   175               case Heap::IN_HEAP:
   157 	for (int i = 0; i < data.size()-1; ++i)
   176                 nd = d + cost[e] + potential[v];
   158 	  (*pot)[data[i].first] += val - data[i].second;
   177                 if (nd < heap[v]) {
   159       }
   178                   heap.decrease(v, nd);
   160 
   179                   pred[v] = e;
   161     }; //class PotentialUpdateMap
   180                 }
   162 
   181                 break;
   163 #ifdef WITH_SCALING
   182               case Heap::POST_HEAP:
   164     /// \brief Map adaptor class for identifing deficit nodes.
   183                 break;
   165     class DeficitBoolMap : public MapBase<ResNode, bool>
   184               }
   166     {
   185             }
   167     private:
   186           }
   168 
   187 
   169       const SupplyRefMap &imb;
   188           // Traversing incoming edges
   170       const Capacity &delta;
   189           for (InEdgeIt e(graph, u); e != INVALID; ++e) {
   171 
   190             if (flow[e] >= delta) {
   172     public:
   191               v = graph.source(e);
   173 
   192               switch(heap.state(v)) {
   174       DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :
   193               case Heap::PRE_HEAP:
   175 	imb(_imb), delta(_delta) {}
   194                 heap.push(v, d - cost[e] + potential[v]);
   176 
   195                 pred[v] = e;
   177       bool operator[](const ResNode &n) const {
   196                 break;
   178 	return imb[n] <= -delta;
   197               case Heap::IN_HEAP:
   179       }
   198                 nd = d - cost[e] + potential[v];
   180 
   199                 if (nd < heap[v]) {
   181     }; //class DeficitBoolMap
   200                   heap.decrease(v, nd);
   182 
   201                   pred[v] = e;
   183     /// \brief Map adaptor class for filtering edges with at least
   202                 }
   184     /// \c delta residual capacity.
   203                 break;
   185     class DeltaFilterMap : public MapBase<ResEdge, bool>
   204               case Heap::POST_HEAP:
   186     {
   205                 break;
   187     private:
   206               }
   188 
   207             }
   189       const ResGraph &gr;
   208           }
   190       const Capacity &delta;
   209         }
   191 
   210         if (heap.empty()) return INVALID;
   192     public:
   211 
   193 
   212         // Updating potentials of processed nodes
   194       DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :
   213         Node t = heap.top();
   195 	gr(_gr), delta(_delta) {}
   214         Cost dt = heap.prio();
   196 
   215         for (int i = 0; i < proc_nodes.size(); ++i)
   197       bool operator[](const ResEdge &e) const {
   216           potential[proc_nodes[i]] -= dist[proc_nodes[i]] - dt;
   198 	return gr.rescap(e) >= delta;
   217 
   199       }
   218         return t;
   200 
   219       }
   201     }; //class DeltaFilterMap
   220 
   202 
   221     }; //class ResidualDijkstra
   203     typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>
       
   204       DeltaResGraph;
       
   205 
       
   206     /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
       
   207     class ResDijkstraTraits :
       
   208       public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>
       
   209     {
       
   210     public:
       
   211 
       
   212       typedef PotentialUpdateMap DistMap;
       
   213 
       
   214       static DistMap *createDistMap(const DeltaResGraph&) {
       
   215 	return new DistMap();
       
   216       }
       
   217 
       
   218     }; //class ResDijkstraTraits
       
   219 
       
   220 #else //WITHOUT_CAPACITY_SCALING
       
   221     /// \brief Map adaptor class for identifing deficit nodes.
       
   222     class DeficitBoolMap : public MapBase<ResNode, bool>
       
   223     {
       
   224     private:
       
   225 
       
   226       const SupplyRefMap &imb;
       
   227 
       
   228     public:
       
   229 
       
   230       DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}
       
   231 
       
   232       bool operator[](const ResNode &n) const {
       
   233 	return imb[n] < 0;
       
   234       }
       
   235 
       
   236     }; //class DeficitBoolMap
       
   237 
       
   238     /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
       
   239     class ResDijkstraTraits :
       
   240       public DijkstraDefaultTraits<ResGraph, ReducedCostMap>
       
   241     {
       
   242     public:
       
   243 
       
   244       typedef PotentialUpdateMap DistMap;
       
   245 
       
   246       static DistMap *createDistMap(const ResGraph&) {
       
   247 	return new DistMap();
       
   248       }
       
   249 
       
   250     }; //class ResDijkstraTraits
       
   251 #endif
       
   252 
   222 
   253   protected:
   223   protected:
   254 
   224 
   255     /// \brief The directed graph the algorithm runs on.
   225     /// \brief The directed graph the algorithm runs on.
   256     const Graph &graph;
   226     const Graph &graph;
   267 
   237 
   268     /// \brief The edge map of the current flow.
   238     /// \brief The edge map of the current flow.
   269     FlowMap flow;
   239     FlowMap flow;
   270     /// \brief The potential node map.
   240     /// \brief The potential node map.
   271     PotentialMap potential;
   241     PotentialMap potential;
   272     /// \brief The residual graph.
   242 
   273     ResGraph res_graph;
   243     /// \brief The residual capacity map.
   274     /// \brief The reduced cost map.
   244     CapacityRefMap res_cap;
   275     ReducedCostMap red_cost;
   245     /// \brief The excess map.
   276 
   246     SupplyRefMap excess;
   277     /// \brief The imbalance map.
   247     /// \brief The excess nodes (i.e. nodes with positive excess).
   278     SupplyRefMap imbalance;
   248     std::vector<Node> excess_nodes;
   279     /// \brief The excess nodes.
       
   280     std::vector<ResNode> excess_nodes;
       
   281     /// \brief The index of the next excess node.
   249     /// \brief The index of the next excess node.
   282     int next_node;
   250     int next_node;
   283 
   251 
   284 #ifdef WITH_SCALING
   252     /// \brief The scaling status (enabled or disabled).
   285     typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>
   253     ScalingEnum scaling;
   286       ResDijkstra;
       
   287     /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
       
   288     /// shortest paths in the residual graph with respect to the
       
   289     /// reduced edge costs.
       
   290     ResDijkstra dijkstra;
       
   291 
       
   292     /// \brief The delta parameter used for capacity scaling.
   254     /// \brief The delta parameter used for capacity scaling.
   293     Capacity delta;
   255     Capacity delta;
   294     /// \brief Edge filter map for filtering edges with at least
   256     /// \brief The maximum number of phases.
   295     /// \c delta residual capacity.
   257     Capacity phase_num;
   296     DeltaFilterMap delta_filter;
       
   297     /// \brief The delta residual graph (i.e. the subgraph of the
       
   298     /// residual graph consisting of edges with at least \c delta
       
   299     /// residual capacity).
       
   300     DeltaResGraph dres_graph;
       
   301     /// \brief Map for identifing deficit nodes.
       
   302     DeficitBoolMap delta_deficit;
       
   303 
       
   304     /// \brief The deficit nodes.
   258     /// \brief The deficit nodes.
   305     std::vector<ResNode> deficit_nodes;
   259     std::vector<Node> deficit_nodes;
   306 
   260 
   307 #else //WITHOUT_CAPACITY_SCALING
   261     /// \brief Implementation of the \ref Dijkstra algorithm for
   308     typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
   262     /// finding augmenting shortest paths in the residual network.
   309       ResDijkstra;
   263     ResidualDijkstra dijkstra;
   310     /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
   264     /// \brief The map of predecessors edges.
   311     /// shortest paths in the residual graph with respect to the
   265     PredMap pred;
   312     /// reduced edge costs.
       
   313     ResDijkstra dijkstra;
       
   314     /// \brief Map for identifing deficit nodes.
       
   315     DeficitBoolMap has_deficit;
       
   316 #endif
       
   317 
       
   318     /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.
       
   319     typename ResDijkstra::PredMap pred;
       
   320     /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to
       
   321     /// update node potentials.
       
   322     PotentialUpdateMap updater;
       
   323 
   266 
   324   public :
   267   public :
   325 
   268 
   326     /// \brief General constructor of the class (with lower bounds).
   269     /// \brief General constructor of the class (with lower bounds).
   327     ///
   270     ///
   331     /// \param _lower The lower bounds of the edges.
   274     /// \param _lower The lower bounds of the edges.
   332     /// \param _capacity The capacities (upper bounds) of the edges.
   275     /// \param _capacity The capacities (upper bounds) of the edges.
   333     /// \param _cost The cost (length) values of the edges.
   276     /// \param _cost The cost (length) values of the edges.
   334     /// \param _supply The supply values of the nodes (signed).
   277     /// \param _supply The supply values of the nodes (signed).
   335     CapacityScaling( const Graph &_graph,
   278     CapacityScaling( const Graph &_graph,
   336 		     const LowerMap &_lower,
   279                      const LowerMap &_lower,
   337 		     const CapacityMap &_capacity,
   280                      const CapacityMap &_capacity,
   338 		     const CostMap &_cost,
   281                      const CostMap &_cost,
   339 		     const SupplyMap &_supply ) :
   282                      const SupplyMap &_supply ) :
   340       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   283       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   341       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   284       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   342       res_graph(_graph, capacity, flow),
   285       res_cap(_graph), excess(_graph), pred(_graph),
   343       red_cost(res_graph, cost, potential), imbalance(_graph),
   286       dijkstra(graph, flow, res_cap, cost, excess, potential, pred)
   344 #ifdef WITH_SCALING
       
   345       delta(0), delta_filter(res_graph, delta),
       
   346       dres_graph(res_graph, delta_filter),
       
   347       dijkstra(dres_graph, red_cost), pred(dres_graph),
       
   348       delta_deficit(imbalance, delta)
       
   349 #else //WITHOUT_CAPACITY_SCALING
       
   350       dijkstra(res_graph, red_cost), pred(res_graph),
       
   351       has_deficit(imbalance)
       
   352 #endif
       
   353     {
   287     {
   354       // Removing nonzero lower bounds
   288       // Removing nonzero lower bounds
   355       capacity = subMap(_capacity, _lower);
   289       capacity = subMap(_capacity, _lower);
       
   290       res_cap = capacity;
   356       Supply sum = 0;
   291       Supply sum = 0;
   357       for (NodeIt n(graph); n != INVALID; ++n) {
   292       for (NodeIt n(graph); n != INVALID; ++n) {
   358 	Supply s = _supply[n];
   293         Supply s = _supply[n];
   359 	for (InEdgeIt e(graph, n); e != INVALID; ++e)
   294         for (InEdgeIt e(graph, n); e != INVALID; ++e)
   360 	  s += _lower[e];
   295           s += _lower[e];
   361 	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
   296         for (OutEdgeIt e(graph, n); e != INVALID; ++e)
   362 	  s -= _lower[e];
   297           s -= _lower[e];
   363 	supply[n] = s;
   298         supply[n] = s;
   364 	sum += s;
   299         sum += s;
   365       }
   300       }
   366       valid_supply = sum == 0;
   301       valid_supply = sum == 0;
   367     }
   302     }
   368 
   303 
   369     /// \brief General constructor of the class (without lower bounds).
   304     /// \brief General constructor of the class (without lower bounds).
   373     /// \param _graph The directed graph the algorithm runs on.
   308     /// \param _graph The directed graph the algorithm runs on.
   374     /// \param _capacity The capacities (upper bounds) of the edges.
   309     /// \param _capacity The capacities (upper bounds) of the edges.
   375     /// \param _cost The cost (length) values of the edges.
   310     /// \param _cost The cost (length) values of the edges.
   376     /// \param _supply The supply values of the nodes (signed).
   311     /// \param _supply The supply values of the nodes (signed).
   377     CapacityScaling( const Graph &_graph,
   312     CapacityScaling( const Graph &_graph,
   378 		     const CapacityMap &_capacity,
   313                      const CapacityMap &_capacity,
   379 		     const CostMap &_cost,
   314                      const CostMap &_cost,
   380 		     const SupplyMap &_supply ) :
   315                      const SupplyMap &_supply ) :
   381       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
   316       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
   382       supply(_supply), flow(_graph, 0), potential(_graph, 0),
   317       supply(_supply), flow(_graph, 0), potential(_graph, 0),
   383       res_graph(_graph, capacity, flow),
   318       res_cap(_capacity), excess(_graph), pred(_graph),
   384       red_cost(res_graph, cost, potential), imbalance(_graph),
   319       dijkstra(graph, flow, res_cap, cost, excess, potential)
   385 #ifdef WITH_SCALING
       
   386       delta(0), delta_filter(res_graph, delta),
       
   387       dres_graph(res_graph, delta_filter),
       
   388       dijkstra(dres_graph, red_cost), pred(dres_graph),
       
   389       delta_deficit(imbalance, delta)
       
   390 #else //WITHOUT_CAPACITY_SCALING
       
   391       dijkstra(res_graph, red_cost), pred(res_graph),
       
   392       has_deficit(imbalance)
       
   393 #endif
       
   394     {
   320     {
   395       // Checking the sum of supply values
   321       // Checking the sum of supply values
   396       Supply sum = 0;
   322       Supply sum = 0;
   397       for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
   323       for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
   398       valid_supply = sum == 0;
   324       valid_supply = sum == 0;
   410     /// \param _t The target node.
   336     /// \param _t The target node.
   411     /// \param _flow_value The required amount of flow from node \c _s
   337     /// \param _flow_value The required amount of flow from node \c _s
   412     /// to node \c _t (i.e. the supply of \c _s and the demand of
   338     /// to node \c _t (i.e. the supply of \c _s and the demand of
   413     /// \c _t).
   339     /// \c _t).
   414     CapacityScaling( const Graph &_graph,
   340     CapacityScaling( const Graph &_graph,
   415 		     const LowerMap &_lower,
   341                      const LowerMap &_lower,
   416 		     const CapacityMap &_capacity,
   342                      const CapacityMap &_capacity,
   417 		     const CostMap &_cost,
   343                      const CostMap &_cost,
   418 		     Node _s, Node _t,
   344                      Node _s, Node _t,
   419 		     Supply _flow_value ) :
   345                      Supply _flow_value ) :
   420       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   346       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   421       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   347       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   422       res_graph(_graph, capacity, flow),
   348       res_cap(_graph), excess(_graph), pred(_graph),
   423       red_cost(res_graph, cost, potential), imbalance(_graph),
   349       dijkstra(graph, flow, res_cap, cost, excess, potential)
   424 #ifdef WITH_SCALING
       
   425       delta(0), delta_filter(res_graph, delta),
       
   426       dres_graph(res_graph, delta_filter),
       
   427       dijkstra(dres_graph, red_cost), pred(dres_graph),
       
   428       delta_deficit(imbalance, delta)
       
   429 #else //WITHOUT_CAPACITY_SCALING
       
   430       dijkstra(res_graph, red_cost), pred(res_graph),
       
   431       has_deficit(imbalance)
       
   432 #endif
       
   433     {
   350     {
   434       // Removing nonzero lower bounds
   351       // Removing nonzero lower bounds
   435       capacity = subMap(_capacity, _lower);
   352       capacity = subMap(_capacity, _lower);
       
   353       res_cap = capacity;
   436       for (NodeIt n(graph); n != INVALID; ++n) {
   354       for (NodeIt n(graph); n != INVALID; ++n) {
   437 	Supply s = 0;
   355         Supply s = 0;
   438 	if (n == _s) s =  _flow_value;
   356         if (n == _s) s =  _flow_value;
   439 	if (n == _t) s = -_flow_value;
   357         if (n == _t) s = -_flow_value;
   440 	for (InEdgeIt e(graph, n); e != INVALID; ++e)
   358         for (InEdgeIt e(graph, n); e != INVALID; ++e)
   441 	  s += _lower[e];
   359           s += _lower[e];
   442 	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
   360         for (OutEdgeIt e(graph, n); e != INVALID; ++e)
   443 	  s -= _lower[e];
   361           s -= _lower[e];
   444 	supply[n] = s;
   362         supply[n] = s;
   445       }
   363       }
   446       valid_supply = true;
   364       valid_supply = true;
   447     }
   365     }
   448 
   366 
   449     /// \brief Simple constructor of the class (without lower bounds).
   367     /// \brief Simple constructor of the class (without lower bounds).
   457     /// \param _t The target node.
   375     /// \param _t The target node.
   458     /// \param _flow_value The required amount of flow from node \c _s
   376     /// \param _flow_value The required amount of flow from node \c _s
   459     /// to node \c _t (i.e. the supply of \c _s and the demand of
   377     /// to node \c _t (i.e. the supply of \c _s and the demand of
   460     /// \c _t).
   378     /// \c _t).
   461     CapacityScaling( const Graph &_graph,
   379     CapacityScaling( const Graph &_graph,
   462 		     const CapacityMap &_capacity,
   380                      const CapacityMap &_capacity,
   463 		     const CostMap &_cost,
   381                      const CostMap &_cost,
   464 		     Node _s, Node _t,
   382                      Node _s, Node _t,
   465 		     Supply _flow_value ) :
   383                      Supply _flow_value ) :
   466       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
   384       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
   467       supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
   385       supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
   468       res_graph(_graph, capacity, flow),
   386       res_cap(_capacity), excess(_graph), pred(_graph),
   469       red_cost(res_graph, cost, potential), imbalance(_graph),
   387       dijkstra(graph, flow, res_cap, cost, excess, potential)
   470 #ifdef WITH_SCALING
       
   471       delta(0), delta_filter(res_graph, delta),
       
   472       dres_graph(res_graph, delta_filter),
       
   473       dijkstra(dres_graph, red_cost), pred(dres_graph),
       
   474       delta_deficit(imbalance, delta)
       
   475 #else //WITHOUT_CAPACITY_SCALING
       
   476       dijkstra(res_graph, red_cost), pred(res_graph),
       
   477       has_deficit(imbalance)
       
   478 #endif
       
   479     {
   388     {
   480       supply[_s] =  _flow_value;
   389       supply[_s] =  _flow_value;
   481       supply[_t] = -_flow_value;
   390       supply[_t] = -_flow_value;
   482       valid_supply = true;
   391       valid_supply = true;
   483     }
   392     }
   509     ///
   418     ///
   510     /// \pre \ref run() must be called before using this function.
   419     /// \pre \ref run() must be called before using this function.
   511     Cost totalCost() const {
   420     Cost totalCost() const {
   512       Cost c = 0;
   421       Cost c = 0;
   513       for (EdgeIt e(graph); e != INVALID; ++e)
   422       for (EdgeIt e(graph); e != INVALID; ++e)
   514 	c += flow[e] * cost[e];
   423         c += flow[e] * cost[e];
   515       return c;
   424       return c;
   516     }
   425     }
   517 
   426 
   518     /// \brief Runs the algorithm.
   427     /// \brief Runs the algorithm.
   519     ///
   428     ///
   520     /// Runs the algorithm.
   429     /// Runs the algorithm.
   521     ///
   430     ///
       
   431     /// \param scaling_mode The scaling mode. In case of WITH_SCALING
       
   432     /// capacity scaling is enabled in the algorithm (this is the
       
   433     /// default value) otherwise it is disabled.
       
   434     /// If the maximum edge capacity and/or the amount of total supply
       
   435     /// is small, the algorithm could be faster without scaling.
       
   436     ///
   522     /// \return \c true if a feasible flow can be found.
   437     /// \return \c true if a feasible flow can be found.
   523     bool run() {
   438     bool run(int scaling_mode = WITH_SCALING) {
   524       return init() && start();
   439       return init(scaling_mode) && start();
   525     }
   440     }
   526 
   441 
   527   protected:
   442   protected:
   528 
   443 
   529     /// \brief Initializes the algorithm.
   444     /// \brief Initializes the algorithm.
   530     bool init() {
   445     bool init(int scaling_mode) {
   531       if (!valid_supply) return false;
   446       if (!valid_supply) return false;
   532       imbalance = supply;
   447       excess = supply;
   533 
   448 
   534       // Initalizing Dijkstra class
       
   535       updater.potentialMap(potential);
       
   536       dijkstra.distMap(updater).predMap(pred);
       
   537 
       
   538 #ifdef WITH_SCALING
       
   539       // Initilaizing delta value
   449       // Initilaizing delta value
   540       Supply max_sup = 0, max_dem = 0;
   450       if (scaling_mode == WITH_SCALING) {
   541       for (NodeIt n(graph); n != INVALID; ++n) {
   451         // With scaling
   542 	if (supply[n] >  max_sup) max_sup =  supply[n];
   452         Supply max_sup = 0, max_dem = 0;
   543 	if (supply[n] < -max_dem) max_dem = -supply[n];
   453         for (NodeIt n(graph); n != INVALID; ++n) {
   544       }
   454           if ( supply[n] > max_sup) max_sup =  supply[n];
   545       if (max_dem < max_sup) max_sup = max_dem;
   455           if (-supply[n] > max_dem) max_dem = -supply[n];
   546       for ( delta = 1; SCALING_FACTOR * delta < max_sup;
   456         }
   547 	    delta *= SCALING_FACTOR ) ;
   457         if (max_dem < max_sup) max_sup = max_dem;
   548 #endif
   458         phase_num = 0;
       
   459         for (delta = 1; 2 * delta <= max_sup; delta *= 2)
       
   460           ++phase_num;
       
   461       } else {
       
   462         // Without scaling
       
   463         delta = 1;
       
   464       }
   549       return true;
   465       return true;
   550     }
   466     }
   551 
   467 
   552 #ifdef WITH_SCALING
   468     /// \brief Executes the algorithm.
       
   469     bool start() {
       
   470       if (delta > 1)
       
   471         return startWithScaling();
       
   472       else
       
   473         return startWithoutScaling();
       
   474     }
       
   475 
   553     /// \brief Executes the capacity scaling version of the successive
   476     /// \brief Executes the capacity scaling version of the successive
   554     /// shortest path algorithm.
   477     /// shortest path algorithm.
   555     bool start() {
   478     bool startWithScaling() {
   556       typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
       
   557       typedef typename DeltaResGraph::Edge DeltaResEdge;
       
   558 #ifdef _DEBUG_ITER_
       
   559       int dijk_num = 0;
       
   560 #endif
       
   561 
       
   562       // Processing capacity scaling phases
   479       // Processing capacity scaling phases
   563       ResNode s, t;
   480       Node s, t;
   564       for ( ; delta >= 1; delta = delta < SCALING_FACTOR && delta > 1 ?
   481       int phase_cnt = 0;
   565 				  1 : delta / SCALING_FACTOR )
   482       int factor = 4;
   566       {
   483       while (true) {
   567 	// Saturating edges not satisfying the optimality condition
   484         // Saturating all edges not satisfying the optimality condition
   568 	Capacity r;
   485         for (EdgeIt e(graph); e != INVALID; ++e) {
   569 	for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {
   486           Node u = graph.source(e), v = graph.target(e);
   570 	  if (red_cost[e] < 0) {
   487           Cost c = cost[e] - potential[u] + potential[v];
   571 	    res_graph.augment(e, r = res_graph.rescap(e));
   488           if (c < 0 && res_cap[e] >= delta) {
   572 	    imbalance[dres_graph.source(e)] -= r;
   489             excess[u] -= res_cap[e];
   573 	    imbalance[dres_graph.target(e)] += r;
   490             excess[v] += res_cap[e];
   574 	  }
   491             flow[e] = capacity[e];
   575 	}
   492             res_cap[e] = 0;
   576 
   493           }
   577 	// Finding excess nodes
   494           else if (c > 0 && flow[e] >= delta) {
   578 	excess_nodes.clear();
   495             excess[u] += flow[e];
   579 	deficit_nodes.clear();
   496             excess[v] -= flow[e];
   580 	for (ResNodeIt n(res_graph); n != INVALID; ++n) {
   497             flow[e] = 0;
   581 	  if (imbalance[n] >=  delta) excess_nodes.push_back(n);
   498             res_cap[e] = capacity[e];
   582 	  if (imbalance[n] <= -delta) deficit_nodes.push_back(n);
   499           }
   583 	}
   500         }
   584 	next_node = 0;
   501 
   585 
   502         // Finding excess nodes and deficit nodes
   586 	// Finding shortest paths
   503         excess_nodes.clear();
   587 	while (next_node < excess_nodes.size()) {
   504         deficit_nodes.clear();
   588 	  // Checking deficit nodes
   505         for (NodeIt n(graph); n != INVALID; ++n) {
   589 	  if (delta > 1) {
   506           if (excess[n] >=  delta) excess_nodes.push_back(n);
   590 	    bool delta_def = false;
   507           if (excess[n] <= -delta) deficit_nodes.push_back(n);
   591 	    for (int i = 0; i < deficit_nodes.size(); ++i) {
   508         }
   592 	      if (imbalance[deficit_nodes[i]] <= -delta) {
   509         next_node = 0;
   593 		delta_def = true;
   510 
   594 		break;
   511         // Finding augmenting shortest paths
   595 	      }
   512         while (next_node < excess_nodes.size()) {
   596 	    }
   513           // Checking deficit nodes
   597 	    if (!delta_def) break;
   514           if (delta > 1) {
   598 	  }
   515             bool delta_deficit = false;
   599 
   516             for (int i = 0; i < deficit_nodes.size(); ++i) {
   600 	  // Running Dijkstra
   517               if (excess[deficit_nodes[i]] <= -delta) {
   601 	  s = excess_nodes[next_node];
   518                 delta_deficit = true;
   602 	  updater.init();
   519                 break;
   603 	  dijkstra.init();
   520               }
   604 	  dijkstra.addSource(s);
   521             }
   605 #ifdef _DEBUG_ITER_
   522             if (!delta_deficit) break;
   606 	  ++dijk_num;
   523           }
   607 #endif
   524 
   608 	  if ((t = dijkstra.start(delta_deficit)) == INVALID) {
   525           // Running Dijkstra
   609 	    if (delta > 1) {
   526           s = excess_nodes[next_node];
   610 	      ++next_node;
   527           if ((t = dijkstra.run(s, delta)) == INVALID) {
   611 	      continue;
   528             if (delta > 1) {
   612 	    }
   529               ++next_node;
   613 	    return false;
   530               continue;
   614 	  }
   531             }
   615 
   532             return false;
   616 	  // Updating node potentials
   533           }
   617 	  updater.update();
   534 
   618 
   535           // Augmenting along a shortest path from s to t.
   619 	  // Augment along a shortest path from s to t
   536           Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
   620 	  Capacity d = imbalance[s] < -imbalance[t] ?
   537           Node u = t;
   621 	    imbalance[s] : -imbalance[t];
   538           Edge e;
   622 	  ResNode u = t;
   539           if (d > delta) {
   623 	  ResEdge e;
   540             while ((e = pred[u]) != INVALID) {
   624 	  if (d > delta) {
   541               Capacity rc;
   625 	    while ((e = pred[u]) != INVALID) {
   542               if (u == graph.target(e)) {
   626 	      if (res_graph.rescap(e) < d) d = res_graph.rescap(e);
   543                 rc = res_cap[e];
   627 	      u = dres_graph.source(e);
   544                 u = graph.source(e);
   628 	    }
   545               } else {
   629 	  }
   546                 rc = flow[e];
   630 	  u = t;
   547                 u = graph.target(e);
   631 	  while ((e = pred[u]) != INVALID) {
   548               }
   632 	    res_graph.augment(e, d);
   549               if (rc < d) d = rc;
   633 	    u = dres_graph.source(e);
   550             }
   634 	  }
   551           }
   635 	  imbalance[s] -= d;
   552           u = t;
   636 	  imbalance[t] += d;
   553           while ((e = pred[u]) != INVALID) {
   637 	  if (imbalance[s] < delta) ++next_node;
   554             if (u == graph.target(e)) {
   638 	}
   555               flow[e] += d;
   639       }
   556               res_cap[e] -= d;
   640 #ifdef _DEBUG_ITER_
   557               u = graph.source(e);
   641       std::cout << "Capacity Scaling algorithm finished with running Dijkstra algorithm "
   558             } else {
   642 	<< dijk_num << " times." << std::endl;
   559               flow[e] -= d;
   643 #endif
   560               res_cap[e] += d;
       
   561               u = graph.target(e);
       
   562             }
       
   563           }
       
   564           excess[s] -= d;
       
   565           excess[t] += d;
       
   566 
       
   567           if (excess[s] < delta) ++next_node;
       
   568         }
       
   569 
       
   570         if (delta == 1) break;
       
   571         if (++phase_cnt > phase_num / 4) factor = 2;
       
   572         delta = delta <= factor ? 1 : delta / factor;
       
   573       }
   644 
   574 
   645       // Handling nonzero lower bounds
   575       // Handling nonzero lower bounds
   646       if (lower) {
   576       if (lower) {
   647 	for (EdgeIt e(graph); e != INVALID; ++e)
   577         for (EdgeIt e(graph); e != INVALID; ++e)
   648 	  flow[e] += (*lower)[e];
   578           flow[e] += (*lower)[e];
   649       }
   579       }
   650       return true;
   580       return true;
   651     }
   581     }
   652 
   582 
   653 #else //WITHOUT_CAPACITY_SCALING
       
   654     /// \brief Executes the successive shortest path algorithm without
   583     /// \brief Executes the successive shortest path algorithm without
   655     /// capacity scaling.
   584     /// capacity scaling.
   656     bool start() {
   585     bool startWithoutScaling() {
   657       // Finding excess nodes
   586       // Finding excess nodes
   658       for (ResNodeIt n(res_graph); n != INVALID; ++n) {
   587       for (NodeIt n(graph); n != INVALID; ++n) {
   659 	if (imbalance[n] > 0) excess_nodes.push_back(n);
   588         if (excess[n] > 0) excess_nodes.push_back(n);
   660       }
   589       }
   661       if (excess_nodes.size() == 0) return true;
   590       if (excess_nodes.size() == 0) return true;
   662       next_node = 0;
   591       next_node = 0;
   663 
   592 
   664       // Finding shortest paths
   593       // Finding shortest paths
   665       ResNode s, t;
   594       Node s, t;
   666       while ( imbalance[excess_nodes[next_node]] > 0 ||
   595       while ( excess[excess_nodes[next_node]] > 0 ||
   667 	      ++next_node < excess_nodes.size() )
   596               ++next_node < excess_nodes.size() )
   668       {
   597       {
   669 	// Running Dijkstra
   598         // Running Dijkstra
   670 	s = excess_nodes[next_node];
   599         s = excess_nodes[next_node];
   671 	updater.init();
   600         if ((t = dijkstra.run(s, 1)) == INVALID)
   672 	dijkstra.init();
   601           return false;
   673 	dijkstra.addSource(s);
   602 
   674 	if ((t = dijkstra.start(has_deficit)) == INVALID)
   603         // Augmenting along a shortest path from s to t
   675 	  return false;
   604         Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
   676 
   605         Node u = t;
   677 	// Updating node potentials
   606         Edge e;
   678 	updater.update();
   607         while ((e = pred[u]) != INVALID) {
   679 
   608           Capacity rc;
   680 	// Augmenting along a shortest path from s to t
   609           if (u == graph.target(e)) {
   681 	Capacity delta = imbalance[s] < -imbalance[t] ?
   610             rc = res_cap[e];
   682 	  imbalance[s] : -imbalance[t];
   611             u = graph.source(e);
   683 	ResNode u = t;
   612           } else {
   684 	ResEdge e;
   613             rc = flow[e];
   685 	while ((e = pred[u]) != INVALID) {
   614             u = graph.target(e);
   686 	  if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);
   615           }
   687 	  u = res_graph.source(e);
   616           if (rc < d) d = rc;
   688 	}
   617         }
   689 	u = t;
   618         u = t;
   690 	while ((e = pred[u]) != INVALID) {
   619         while ((e = pred[u]) != INVALID) {
   691 	  res_graph.augment(e, delta);
   620           if (u == graph.target(e)) {
   692 	  u = res_graph.source(e);
   621             flow[e] += d;
   693 	}
   622             res_cap[e] -= d;
   694 	imbalance[s] -= delta;
   623             u = graph.source(e);
   695 	imbalance[t] += delta;
   624           } else {
       
   625             flow[e] -= d;
       
   626             res_cap[e] += d;
       
   627             u = graph.target(e);
       
   628           }
       
   629         }
       
   630         excess[s] -= d;
       
   631         excess[t] += d;
   696       }
   632       }
   697 
   633 
   698       // Handling nonzero lower bounds
   634       // Handling nonzero lower bounds
   699       if (lower) {
   635       if (lower) {
   700 	for (EdgeIt e(graph); e != INVALID; ++e)
   636         for (EdgeIt e(graph); e != INVALID; ++e)
   701 	  flow[e] += (*lower)[e];
   637           flow[e] += (*lower)[e];
   702       }
   638       }
   703       return true;
   639       return true;
   704     }
   640     }
   705 #endif
       
   706 
   641 
   707   }; //class CapacityScaling
   642   }; //class CapacityScaling
   708 
   643 
   709   ///@}
   644   ///@}
   710 
   645