lemon/capacity_scaling.h
changeset 2557 673cb4d1060b
parent 2553 bfced05fa852
child 2574 7058c9690e7d
equal deleted inserted replaced
7:a9c67524c45e 8:985ac162a9a2
    20 #define LEMON_CAPACITY_SCALING_H
    20 #define LEMON_CAPACITY_SCALING_H
    21 
    21 
    22 /// \ingroup min_cost_flow
    22 /// \ingroup min_cost_flow
    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 flow.
    26 /// flow.
       
    27 
    26 
    28 #include <lemon/graph_adaptor.h>
    27 #include <lemon/graph_adaptor.h>
    29 #include <lemon/bin_heap.h>
    28 #include <lemon/bin_heap.h>
    30 #include <vector>
    29 #include <vector>
    31 
    30 
    47   /// \param CapacityMap The type of the capacity (upper bound) map.
    46   /// \param CapacityMap The type of the capacity (upper bound) map.
    48   /// \param CostMap The type of the cost (length) map.
    47   /// \param CostMap The type of the cost (length) map.
    49   /// \param SupplyMap The type of the supply map.
    48   /// \param SupplyMap The type of the supply map.
    50   ///
    49   ///
    51   /// \warning
    50   /// \warning
    52   /// - Edge capacities and costs should be nonnegative integers.
    51   /// - Edge capacities and costs should be non-negative integers.
    53   ///   However \c CostMap::Value should be signed type.
    52   ///   However \c CostMap::Value should be signed type.
    54   /// - Supply values should be signed integers.
    53   /// - Supply values should be signed integers.
    55   /// - \c LowerMap::Value must be convertible to
    54   /// - \c LowerMap::Value must be convertible to
    56   ///   \c CapacityMap::Value and \c CapacityMap::Value must be
    55   ///   \c CapacityMap::Value and \c CapacityMap::Value must be
    57   ///   convertible to \c SupplyMap::Value.
    56   ///   convertible to \c SupplyMap::Value.
    64              typename CostMap = typename Graph::template EdgeMap<int>,
    63              typename CostMap = typename Graph::template EdgeMap<int>,
    65              typename SupplyMap = typename Graph::template NodeMap
    64              typename SupplyMap = typename Graph::template NodeMap
    66                                   <typename CapacityMap::Value> >
    65                                   <typename CapacityMap::Value> >
    67   class CapacityScaling
    66   class CapacityScaling
    68   {
    67   {
    69     typedef typename Graph::Node Node;
    68     GRAPH_TYPEDEFS(typename Graph);
    70     typedef typename Graph::NodeIt NodeIt;
       
    71     typedef typename Graph::Edge Edge;
       
    72     typedef typename Graph::EdgeIt EdgeIt;
       
    73     typedef typename Graph::InEdgeIt InEdgeIt;
       
    74     typedef typename Graph::OutEdgeIt OutEdgeIt;
       
    75 
    69 
    76     typedef typename LowerMap::Value Lower;
    70     typedef typename LowerMap::Value Lower;
    77     typedef typename CapacityMap::Value Capacity;
    71     typedef typename CapacityMap::Value Capacity;
    78     typedef typename CostMap::Value Cost;
    72     typedef typename CostMap::Value Cost;
    79     typedef typename SupplyMap::Value Supply;
    73     typedef typename SupplyMap::Value Supply;
    80     typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
    74     typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
    81     typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
    75     typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
    82     typedef typename Graph::template NodeMap<Edge> PredMap;
    76     typedef typename Graph::template NodeMap<Edge> PredMap;
    83 
    77 
    84   public:
    78   public:
    85 
    79 
    86     /// \brief Type to enable or disable capacity scaling.
    80     /// Type to enable or disable capacity scaling.
    87     enum ScalingEnum {
    81     enum ScalingEnum {
    88       WITH_SCALING = 0,
    82       WITH_SCALING = 0,
    89       WITHOUT_SCALING = -1
    83       WITHOUT_SCALING = -1
    90     };
    84     };
    91 
    85 
    92     /// \brief The type of the flow map.
    86     /// The type of the flow map.
    93     typedef CapacityRefMap FlowMap;
    87     typedef typename Graph::template EdgeMap<Capacity> FlowMap;
    94     /// \brief The type of the potential map.
    88     /// The type of the potential map.
    95     typedef typename Graph::template NodeMap<Cost> PotentialMap;
    89     typedef typename Graph::template NodeMap<Cost> PotentialMap;
    96 
    90 
    97   protected:
    91   protected:
    98 
    92 
    99     /// \brief Special implementation of the \ref Dijkstra algorithm
    93     /// \brief Special implementation of the \ref Dijkstra algorithm
   108       typedef typename Graph::template NodeMap<int> HeapCrossRef;
   102       typedef typename Graph::template NodeMap<int> HeapCrossRef;
   109       typedef BinHeap<Cost, HeapCrossRef> Heap;
   103       typedef BinHeap<Cost, HeapCrossRef> Heap;
   110 
   104 
   111     protected:
   105     protected:
   112 
   106 
   113       /// \brief The directed graph the algorithm runs on.
   107       /// The directed graph the algorithm runs on.
   114       const Graph &graph;
   108       const Graph &graph;
   115 
   109 
   116       /// \brief The flow map.
   110       /// The flow map.
   117       const FlowMap &flow;
   111       const FlowMap &flow;
   118       /// \brief The residual capacity map.
   112       /// The residual capacity map.
   119       const CapacityRefMap &res_cap;
   113       const CapacityEdgeMap &res_cap;
   120       /// \brief The cost map.
   114       /// The cost map.
   121       const CostMap &cost;
   115       const CostMap &cost;
   122       /// \brief The excess map.
   116       /// The excess map.
   123       const SupplyRefMap &excess;
   117       const SupplyNodeMap &excess;
   124 
   118 
   125       /// \brief The potential map.
   119       /// The potential map.
   126       PotentialMap &potential;
   120       PotentialMap &potential;
   127 
   121 
   128       /// \brief The distance map.
   122       /// The distance map.
   129       CostNodeMap dist;
   123       CostNodeMap dist;
   130       /// \brief The map of predecessors edges.
   124       /// The map of predecessors edges.
   131       PredMap &pred;
   125       PredMap &pred;
   132       /// \brief The processed (i.e. permanently labeled) nodes.
   126       /// The processed (i.e. permanently labeled) nodes.
   133       std::vector<Node> proc_nodes;
   127       std::vector<Node> proc_nodes;
   134 
   128 
   135     public:
   129     public:
   136 
   130 
   137       /// \brief The constructor of the class.
   131       /// The constructor of the class.
   138       ResidualDijkstra( const Graph &_graph,
   132       ResidualDijkstra( const Graph &_graph,
   139                         const FlowMap &_flow,
   133                         const FlowMap &_flow,
   140                         const CapacityRefMap &_res_cap,
   134                         const CapacityEdgeMap &_res_cap,
   141                         const CostMap &_cost,
   135                         const CostMap &_cost,
   142                         const SupplyMap &_excess,
   136                         const SupplyMap &_excess,
   143                         PotentialMap &_potential,
   137                         PotentialMap &_potential,
   144                         PredMap &_pred ) :
   138                         PredMap &_pred ) :
   145         graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost),
   139         graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost),
   146         excess(_excess), potential(_potential), dist(_graph),
   140         excess(_excess), potential(_potential), dist(_graph),
   147         pred(_pred)
   141         pred(_pred)
   148       {}
   142       {}
   149 
   143 
   150       /// \brief Runs the algorithm from the given source node.
   144       /// Runs the algorithm from the given source node.
   151       Node run(Node s, Capacity delta) {
   145       Node run(Node s, Capacity delta) {
   152         HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP);
   146         HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP);
   153         Heap heap(heap_cross_ref);
   147         Heap heap(heap_cross_ref);
   154         heap.push(s, 0);
   148         heap.push(s, 0);
   155         pred[s] = INVALID;
   149         pred[s] = INVALID;
   220 
   214 
   221     }; //class ResidualDijkstra
   215     }; //class ResidualDijkstra
   222 
   216 
   223   protected:
   217   protected:
   224 
   218 
   225     /// \brief The directed graph the algorithm runs on.
   219     /// The directed graph the algorithm runs on.
   226     const Graph &graph;
   220     const Graph &graph;
   227     /// \brief The original lower bound map.
   221     /// The original lower bound map.
   228     const LowerMap *lower;
   222     const LowerMap *lower;
   229     /// \brief The modified capacity map.
   223     /// The modified capacity map.
   230     CapacityRefMap capacity;
   224     CapacityEdgeMap capacity;
   231     /// \brief The cost map.
   225     /// The cost map.
   232     const CostMap &cost;
   226     const CostMap &cost;
   233     /// \brief The modified supply map.
   227     /// The modified supply map.
   234     SupplyRefMap supply;
   228     SupplyNodeMap supply;
   235     /// \brief The sum of supply values equals zero.
       
   236     bool valid_supply;
   229     bool valid_supply;
   237 
   230 
   238     /// \brief The edge map of the current flow.
   231     /// The edge map of the current flow.
   239     FlowMap flow;
   232     FlowMap flow;
   240     /// \brief The potential node map.
   233     /// The potential node map.
   241     PotentialMap potential;
   234     PotentialMap potential;
   242 
   235 
   243     /// \brief The residual capacity map.
   236     /// The residual capacity map.
   244     CapacityRefMap res_cap;
   237     CapacityEdgeMap res_cap;
   245     /// \brief The excess map.
   238     /// The excess map.
   246     SupplyRefMap excess;
   239     SupplyNodeMap excess;
   247     /// \brief The excess nodes (i.e. nodes with positive excess).
   240     /// The excess nodes (i.e. the nodes with positive excess).
   248     std::vector<Node> excess_nodes;
   241     std::vector<Node> excess_nodes;
   249     /// \brief The index of the next excess node.
   242     /// The deficit nodes (i.e. the nodes with negative excess).
   250     int next_node;
   243     std::vector<Node> deficit_nodes;
   251 
   244 
   252     /// \brief The scaling status (enabled or disabled).
   245     /// The scaling status (enabled or disabled).
   253     ScalingEnum scaling;
   246     ScalingEnum scaling;
   254     /// \brief The delta parameter used for capacity scaling.
   247     /// The \c delta parameter used for capacity scaling.
   255     Capacity delta;
   248     Capacity delta;
   256     /// \brief The maximum number of phases.
   249     /// The maximum number of phases.
   257     Capacity phase_num;
   250     int phase_num;
   258     /// \brief The deficit nodes.
       
   259     std::vector<Node> deficit_nodes;
       
   260 
   251 
   261     /// \brief Implementation of the \ref Dijkstra algorithm for
   252     /// \brief Implementation of the \ref Dijkstra algorithm for
   262     /// finding augmenting shortest paths in the residual network.
   253     /// finding augmenting shortest paths in the residual network.
   263     ResidualDijkstra dijkstra;
   254     ResidualDijkstra dijkstra;
   264     /// \brief The map of predecessors edges.
   255     /// The map of predecessors edges.
   265     PredMap pred;
   256     PredMap pred;
   266 
   257 
   267   public :
   258   public :
   268 
   259 
   269     /// \brief General constructor of the class (with lower bounds).
   260     /// \brief General constructor of the class (with lower bounds).
   283       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   274       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   284       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   275       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   285       res_cap(_graph), excess(_graph), pred(_graph),
   276       res_cap(_graph), excess(_graph), pred(_graph),
   286       dijkstra(graph, flow, res_cap, cost, excess, potential, pred)
   277       dijkstra(graph, flow, res_cap, cost, excess, potential, pred)
   287     {
   278     {
   288       // Removing nonzero lower bounds
   279       // Removing non-zero lower bounds
   289       capacity = subMap(_capacity, _lower);
   280       capacity = subMap(_capacity, _lower);
   290       res_cap = capacity;
   281       res_cap = capacity;
   291       Supply sum = 0;
   282       Supply sum = 0;
   292       for (NodeIt n(graph); n != INVALID; ++n) {
   283       for (NodeIt n(graph); n != INVALID; ++n) {
   293         Supply s = _supply[n];
   284         Supply s = _supply[n];
   333     /// \param _capacity The capacities (upper bounds) of the edges.
   324     /// \param _capacity The capacities (upper bounds) of the edges.
   334     /// \param _cost The cost (length) values of the edges.
   325     /// \param _cost The cost (length) values of the edges.
   335     /// \param _s The source node.
   326     /// \param _s The source node.
   336     /// \param _t The target node.
   327     /// \param _t The target node.
   337     /// \param _flow_value The required amount of flow from node \c _s
   328     /// \param _flow_value The required amount of flow from node \c _s
   338     /// to node \c _t (i.e. the supply of \c _s and the demand of
   329     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
   339     /// \c _t).
       
   340     CapacityScaling( const Graph &_graph,
   330     CapacityScaling( const Graph &_graph,
   341                      const LowerMap &_lower,
   331                      const LowerMap &_lower,
   342                      const CapacityMap &_capacity,
   332                      const CapacityMap &_capacity,
   343                      const CostMap &_cost,
   333                      const CostMap &_cost,
   344                      Node _s, Node _t,
   334                      Node _s, Node _t,
   346       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   336       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
   347       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   337       supply(_graph), flow(_graph, 0), potential(_graph, 0),
   348       res_cap(_graph), excess(_graph), pred(_graph),
   338       res_cap(_graph), excess(_graph), pred(_graph),
   349       dijkstra(graph, flow, res_cap, cost, excess, potential)
   339       dijkstra(graph, flow, res_cap, cost, excess, potential)
   350     {
   340     {
   351       // Removing nonzero lower bounds
   341       // Removing non-zero lower bounds
   352       capacity = subMap(_capacity, _lower);
   342       capacity = subMap(_capacity, _lower);
   353       res_cap = capacity;
   343       res_cap = capacity;
   354       for (NodeIt n(graph); n != INVALID; ++n) {
   344       for (NodeIt n(graph); n != INVALID; ++n) {
   355         Supply s = 0;
   345         Supply s = 0;
   356         if (n == _s) s =  _flow_value;
   346         if (n == _s) s =  _flow_value;
   372     /// \param _capacity The capacities (upper bounds) of the edges.
   362     /// \param _capacity The capacities (upper bounds) of the edges.
   373     /// \param _cost The cost (length) values of the edges.
   363     /// \param _cost The cost (length) values of the edges.
   374     /// \param _s The source node.
   364     /// \param _s The source node.
   375     /// \param _t The target node.
   365     /// \param _t The target node.
   376     /// \param _flow_value The required amount of flow from node \c _s
   366     /// \param _flow_value The required amount of flow from node \c _s
   377     /// to node \c _t (i.e. the supply of \c _s and the demand of
   367     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
   378     /// \c _t).
       
   379     CapacityScaling( const Graph &_graph,
   368     CapacityScaling( const Graph &_graph,
   380                      const CapacityMap &_capacity,
   369                      const CapacityMap &_capacity,
   381                      const CostMap &_cost,
   370                      const CostMap &_cost,
   382                      Node _s, Node _t,
   371                      Node _s, Node _t,
   383                      Supply _flow_value ) :
   372                      Supply _flow_value ) :
   389       supply[_s] =  _flow_value;
   378       supply[_s] =  _flow_value;
   390       supply[_t] = -_flow_value;
   379       supply[_t] = -_flow_value;
   391       valid_supply = true;
   380       valid_supply = true;
   392     }
   381     }
   393 
   382 
       
   383     /// \brief Runs the algorithm.
       
   384     ///
       
   385     /// Runs the algorithm.
       
   386     ///
       
   387     /// \param scaling_mode The scaling mode. In case of WITH_SCALING
       
   388     /// capacity scaling is enabled in the algorithm (this is the
       
   389     /// default) otherwise it is disabled.
       
   390     /// If the maximum edge capacity and/or the amount of total supply
       
   391     /// is small, the algorithm could be slightly faster without
       
   392     /// scaling.
       
   393     ///
       
   394     /// \return \c true if a feasible flow can be found.
       
   395     bool run(int scaling_mode = WITH_SCALING) {
       
   396       return init(scaling_mode) && start();
       
   397     }
       
   398 
   394     /// \brief Returns a const reference to the flow map.
   399     /// \brief Returns a const reference to the flow map.
   395     ///
   400     ///
   396     /// Returns a const reference to the flow map.
   401     /// Returns a const reference to the flow map.
   397     ///
   402     ///
   398     /// \pre \ref run() must be called before using this function.
   403     /// \pre \ref run() must be called before using this function.
   422       for (EdgeIt e(graph); e != INVALID; ++e)
   427       for (EdgeIt e(graph); e != INVALID; ++e)
   423         c += flow[e] * cost[e];
   428         c += flow[e] * cost[e];
   424       return c;
   429       return c;
   425     }
   430     }
   426 
   431 
   427     /// \brief Runs the algorithm.
       
   428     ///
       
   429     /// Runs the algorithm.
       
   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     ///
       
   437     /// \return \c true if a feasible flow can be found.
       
   438     bool run(int scaling_mode = WITH_SCALING) {
       
   439       return init(scaling_mode) && start();
       
   440     }
       
   441 
       
   442   protected:
   432   protected:
   443 
   433 
   444     /// \brief Initializes the algorithm.
   434     /// Initializes the algorithm.
   445     bool init(int scaling_mode) {
   435     bool init(int scaling_mode) {
   446       if (!valid_supply) return false;
   436       if (!valid_supply) return false;
   447       excess = supply;
   437       excess = supply;
   448 
   438 
   449       // Initilaizing delta value
   439       // Initilaizing delta value
   463         delta = 1;
   453         delta = 1;
   464       }
   454       }
   465       return true;
   455       return true;
   466     }
   456     }
   467 
   457 
   468     /// \brief Executes the algorithm.
   458     /// Executes the algorithm.
   469     bool start() {
   459     bool start() {
   470       if (delta > 1)
   460       if (delta > 1)
   471         return startWithScaling();
   461         return startWithScaling();
   472       else
   462       else
   473         return startWithoutScaling();
   463         return startWithoutScaling();
   504         deficit_nodes.clear();
   494         deficit_nodes.clear();
   505         for (NodeIt n(graph); n != INVALID; ++n) {
   495         for (NodeIt n(graph); n != INVALID; ++n) {
   506           if (excess[n] >=  delta) excess_nodes.push_back(n);
   496           if (excess[n] >=  delta) excess_nodes.push_back(n);
   507           if (excess[n] <= -delta) deficit_nodes.push_back(n);
   497           if (excess[n] <= -delta) deficit_nodes.push_back(n);
   508         }
   498         }
   509         next_node = 0;
   499         int next_node = 0;
   510 
   500 
   511         // Finding augmenting shortest paths
   501         // Finding augmenting shortest paths
   512         while (next_node < excess_nodes.size()) {
   502         while (next_node < excess_nodes.size()) {
   513           // Checking deficit nodes
   503           // Checking deficit nodes
   514           if (delta > 1) {
   504           if (delta > 1) {
   570         if (delta == 1) break;
   560         if (delta == 1) break;
   571         if (++phase_cnt > phase_num / 4) factor = 2;
   561         if (++phase_cnt > phase_num / 4) factor = 2;
   572         delta = delta <= factor ? 1 : delta / factor;
   562         delta = delta <= factor ? 1 : delta / factor;
   573       }
   563       }
   574 
   564 
   575       // Handling nonzero lower bounds
   565       // Handling non-zero lower bounds
   576       if (lower) {
   566       if (lower) {
   577         for (EdgeIt e(graph); e != INVALID; ++e)
   567         for (EdgeIt e(graph); e != INVALID; ++e)
   578           flow[e] += (*lower)[e];
   568           flow[e] += (*lower)[e];
   579       }
   569       }
   580       return true;
   570       return true;
   582 
   572 
   583     /// \brief Executes the successive shortest path algorithm without
   573     /// \brief Executes the successive shortest path algorithm without
   584     /// capacity scaling.
   574     /// capacity scaling.
   585     bool startWithoutScaling() {
   575     bool startWithoutScaling() {
   586       // Finding excess nodes
   576       // Finding excess nodes
   587       for (NodeIt n(graph); n != INVALID; ++n) {
   577       for (NodeIt n(graph); n != INVALID; ++n)
   588         if (excess[n] > 0) excess_nodes.push_back(n);
   578         if (excess[n] > 0) excess_nodes.push_back(n);
   589       }
       
   590       if (excess_nodes.size() == 0) return true;
   579       if (excess_nodes.size() == 0) return true;
   591       next_node = 0;
   580       int next_node = 0;
   592 
   581 
   593       // Finding shortest paths
   582       // Finding shortest paths
   594       Node s, t;
   583       Node s, t;
   595       while ( excess[excess_nodes[next_node]] > 0 ||
   584       while ( excess[excess_nodes[next_node]] > 0 ||
   596               ++next_node < excess_nodes.size() )
   585               ++next_node < excess_nodes.size() )
   629         }
   618         }
   630         excess[s] -= d;
   619         excess[s] -= d;
   631         excess[t] += d;
   620         excess[t] += d;
   632       }
   621       }
   633 
   622 
   634       // Handling nonzero lower bounds
   623       // Handling non-zero lower bounds
   635       if (lower) {
   624       if (lower) {
   636         for (EdgeIt e(graph); e != INVALID; ++e)
   625         for (EdgeIt e(graph); e != INVALID; ++e)
   637           flow[e] += (*lower)[e];
   626           flow[e] += (*lower)[e];
   638       }
   627       }
   639       return true;
   628       return true;