lemon/capacity_scaling.h
changeset 2588 4d3bc1d04c1d
parent 2581 054566ac0934
child 2589 1bbb28acb8c9
equal deleted inserted replaced
10:56918ba44652 11:38ed39126590
    23 ///
    23 ///
    24 /// \file
    24 /// \file
    25 /// \brief Capacity scaling algorithm for finding a minimum cost flow.
    25 /// \brief Capacity scaling algorithm for finding a minimum cost flow.
    26 
    26 
    27 #include <vector>
    27 #include <vector>
    28 
       
    29 #include <lemon/graph_adaptor.h>
       
    30 #include <lemon/bin_heap.h>
    28 #include <lemon/bin_heap.h>
    31 
    29 
    32 namespace lemon {
    30 namespace lemon {
    33 
    31 
    34   /// \addtogroup min_cost_flow
    32   /// \addtogroup min_cost_flow
    88     /// residual network of the graph with respect to the reduced edge
    86     /// residual network of the graph with respect to the reduced edge
    89     /// costs and modifying the node potentials according to the
    87     /// costs and modifying the node potentials according to the
    90     /// distance of the nodes.
    88     /// distance of the nodes.
    91     class ResidualDijkstra
    89     class ResidualDijkstra
    92     {
    90     {
    93       typedef typename Graph::template NodeMap<Cost> CostNodeMap;
       
    94       typedef typename Graph::template NodeMap<Edge> PredMap;
       
    95 
       
    96       typedef typename Graph::template NodeMap<int> HeapCrossRef;
    91       typedef typename Graph::template NodeMap<int> HeapCrossRef;
    97       typedef BinHeap<Cost, HeapCrossRef> Heap;
    92       typedef BinHeap<Cost, HeapCrossRef> Heap;
    98 
    93 
    99     private:
    94     private:
   100 
    95 
   107       const CostMap &_cost;
   102       const CostMap &_cost;
   108       const SupplyNodeMap &_excess;
   103       const SupplyNodeMap &_excess;
   109       PotentialMap &_potential;
   104       PotentialMap &_potential;
   110 
   105 
   111       // The distance map
   106       // The distance map
   112       CostNodeMap _dist;
   107       PotentialMap _dist;
   113       // The pred edge map
   108       // The pred edge map
   114       PredMap &_pred;
   109       PredMap &_pred;
   115       // The processed (i.e. permanently labeled) nodes
   110       // The processed (i.e. permanently labeled) nodes
   116       std::vector<Node> _proc_nodes;
   111       std::vector<Node> _proc_nodes;
   117 
   112 
   129         _excess(excess), _potential(potential), _dist(graph),
   124         _excess(excess), _potential(potential), _dist(graph),
   130         _pred(pred)
   125         _pred(pred)
   131       {}
   126       {}
   132 
   127 
   133       /// Runs the algorithm from the given source node.
   128       /// Runs the algorithm from the given source node.
   134       Node run(Node s, Capacity delta) {
   129       Node run(Node s, Capacity delta = 1) {
   135         HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
   130         HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
   136         Heap heap(heap_cross_ref);
   131         Heap heap(heap_cross_ref);
   137         heap.push(s, 0);
   132         heap.push(s, 0);
   138         _pred[s] = INVALID;
   133         _pred[s] = INVALID;
   139         _proc_nodes.clear();
   134         _proc_nodes.clear();
   452     /// \pre \ref run() must be called before using this function.
   447     /// \pre \ref run() must be called before using this function.
   453     const PotentialMap& potentialMap() const {
   448     const PotentialMap& potentialMap() const {
   454       return *_potential;
   449       return *_potential;
   455     }
   450     }
   456 
   451 
   457     /// \brief Returns the flow on the edge.
   452     /// \brief Returns the flow on the given edge.
   458     ///
   453     ///
   459     /// Returns the flow on the edge.
   454     /// Returns the flow on the given edge.
   460     ///
   455     ///
   461     /// \pre \ref run() must be called before using this function.
   456     /// \pre \ref run() must be called before using this function.
   462     Capacity flow(const Edge& edge) const {
   457     Capacity flow(const Edge& edge) const {
   463       return (*_flow)[edge];
   458       return (*_flow)[edge];
   464     }
   459     }
   465 
   460 
   466     /// \brief Returns the potential of the node.
   461     /// \brief Returns the potential of the given node.
   467     ///
   462     ///
   468     /// Returns the potential of the node.
   463     /// Returns the potential of the given node.
   469     ///
   464     ///
   470     /// \pre \ref run() must be called before using this function.
   465     /// \pre \ref run() must be called before using this function.
   471     Cost potential(const Node& node) const {
   466     Cost potential(const Node& node) const {
   472       return (*_potential)[node];
   467       return (*_potential)[node];
   473     }
   468     }
   515         Supply max_sup = 0, max_dem = 0;
   510         Supply max_sup = 0, max_dem = 0;
   516         for (NodeIt n(_graph); n != INVALID; ++n) {
   511         for (NodeIt n(_graph); n != INVALID; ++n) {
   517           if ( _supply[n] > max_sup) max_sup =  _supply[n];
   512           if ( _supply[n] > max_sup) max_sup =  _supply[n];
   518           if (-_supply[n] > max_dem) max_dem = -_supply[n];
   513           if (-_supply[n] > max_dem) max_dem = -_supply[n];
   519         }
   514         }
   520         if (max_dem < max_sup) max_sup = max_dem;
   515         Capacity max_cap = 0;
       
   516         for (EdgeIt e(_graph); e != INVALID; ++e) {
       
   517           if (_capacity[e] > max_cap) max_cap = _capacity[e];
       
   518         }
       
   519         max_sup = std::min(std::min(max_sup, max_dem), max_cap);
   521         _phase_num = 0;
   520         _phase_num = 0;
   522         for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
   521         for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
   523           ++_phase_num;
   522           ++_phase_num;
   524       } else {
   523       } else {
   525         // Without scaling
   524         // Without scaling
   593             }
   592             }
   594             return false;
   593             return false;
   595           }
   594           }
   596 
   595 
   597           // Augmenting along a shortest path from s to t.
   596           // Augmenting along a shortest path from s to t.
   598           Capacity d = _excess[s] < -_excess[t] ? _excess[s] : -_excess[t];
   597           Capacity d = std::min(_excess[s], -_excess[t]);
   599           Node u = t;
   598           Node u = t;
   600           Edge e;
   599           Edge e;
   601           if (d > _delta) {
   600           if (d > _delta) {
   602             while ((e = _pred[u]) != INVALID) {
   601             while ((e = _pred[u]) != INVALID) {
   603               Capacity rc;
   602               Capacity rc;
   655       while ( _excess[_excess_nodes[next_node]] > 0 ||
   654       while ( _excess[_excess_nodes[next_node]] > 0 ||
   656               ++next_node < int(_excess_nodes.size()) )
   655               ++next_node < int(_excess_nodes.size()) )
   657       {
   656       {
   658         // Running Dijkstra
   657         // Running Dijkstra
   659         s = _excess_nodes[next_node];
   658         s = _excess_nodes[next_node];
   660         if ((t = _dijkstra->run(s, 1)) == INVALID)
   659         if ((t = _dijkstra->run(s)) == INVALID) break;
   661           return false;
       
   662 
   660 
   663         // Augmenting along a shortest path from s to t
   661         // Augmenting along a shortest path from s to t
   664         Capacity d = _excess[s] < -_excess[t] ? _excess[s] : -_excess[t];
   662         Capacity d = std::min(_excess[s], -_excess[t]);
   665         Node u = t;
   663         Node u = t;
   666         Edge e;
   664         Edge e;
   667         while ((e = _pred[u]) != INVALID) {
   665         if (d > 1) {
   668           Capacity rc;
   666           while ((e = _pred[u]) != INVALID) {
   669           if (u == _graph.target(e)) {
   667             Capacity rc;
   670             rc = _res_cap[e];
   668             if (u == _graph.target(e)) {
   671             u = _graph.source(e);
   669               rc = _res_cap[e];
   672           } else {
   670               u = _graph.source(e);
   673             rc = (*_flow)[e];
   671             } else {
   674             u = _graph.target(e);
   672               rc = (*_flow)[e];
   675           }
   673               u = _graph.target(e);
   676           if (rc < d) d = rc;
   674             }
       
   675             if (rc < d) d = rc;
       
   676           }
   677         }
   677         }
   678         u = t;
   678         u = t;
   679         while ((e = _pred[u]) != INVALID) {
   679         while ((e = _pred[u]) != INVALID) {
   680           if (u == _graph.target(e)) {
   680           if (u == _graph.target(e)) {
   681             (*_flow)[e] += d;
   681             (*_flow)[e] += d;