COIN-OR::LEMON - Graph Library

Changeset 809:22bb98ca0101 in lemon-1.2 for lemon/cost_scaling.h


Ignore:
Timestamp:
11/12/09 23:30:45 (10 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Rebase:
35333461653237666330353538396562623465656466303531646265343536366333386130366130
Message:

Entirely rework CostScaling? (#180)

  • Use the new interface similarly to NetworkSimplex?.
  • Rework the implementation using an efficient internal structure for handling the residual network. This improvement made the code much faster.
  • Handle GEQ supply type (LEQ is not supported).
  • Handle infinite upper bounds.
  • Handle negative costs (for arcs of finite upper bound).
  • Traits class + named parameter for the LargeCost? type used in internal computations.
  • Extend the documentation.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/cost_scaling.h

    r808 r809  
    3131#include <lemon/maps.h>
    3232#include <lemon/math.h>
    33 #include <lemon/adaptors.h>
     33#include <lemon/static_graph.h>
    3434#include <lemon/circulation.h>
    3535#include <lemon/bellman_ford.h>
     
    3737namespace lemon {
    3838
     39  /// \brief Default traits class of CostScaling algorithm.
     40  ///
     41  /// Default traits class of CostScaling algorithm.
     42  /// \tparam GR Digraph type.
     43  /// \tparam V The value type used for flow amounts, capacity bounds
     44  /// and supply values. By default it is \c int.
     45  /// \tparam C The value type used for costs and potentials.
     46  /// By default it is the same as \c V.
     47#ifdef DOXYGEN
     48  template <typename GR, typename V = int, typename C = V>
     49#else
     50  template < typename GR, typename V = int, typename C = V,
     51             bool integer = std::numeric_limits<C>::is_integer >
     52#endif
     53  struct CostScalingDefaultTraits
     54  {
     55    /// The type of the digraph
     56    typedef GR Digraph;
     57    /// The type of the flow amounts, capacity bounds and supply values
     58    typedef V Value;
     59    /// The type of the arc costs
     60    typedef C Cost;
     61
     62    /// \brief The large cost type used for internal computations
     63    ///
     64    /// The large cost type used for internal computations.
     65    /// It is \c long \c long if the \c Cost type is integer,
     66    /// otherwise it is \c double.
     67    /// \c Cost must be convertible to \c LargeCost.
     68    typedef double LargeCost;
     69  };
     70
     71  // Default traits class for integer cost types
     72  template <typename GR, typename V, typename C>
     73  struct CostScalingDefaultTraits<GR, V, C, true>
     74  {
     75    typedef GR Digraph;
     76    typedef V Value;
     77    typedef C Cost;
     78#ifdef LEMON_HAVE_LONG_LONG
     79    typedef long long LargeCost;
     80#else
     81    typedef long LargeCost;
     82#endif
     83  };
     84
     85
    3986  /// \addtogroup min_cost_flow_algs
    4087  /// @{
    4188
    42   /// \brief Implementation of the cost scaling algorithm for finding a
    43   /// minimum cost flow.
     89  /// \brief Implementation of the Cost Scaling algorithm for
     90  /// finding a \ref min_cost_flow "minimum cost flow".
    4491  ///
    45   /// \ref CostScaling implements the cost scaling algorithm performing
    46   /// augment/push and relabel operations for finding a minimum cost
    47   /// flow.
     92  /// \ref CostScaling implements a cost scaling algorithm that performs
     93  /// push/augment and relabel operations for finding a minimum cost
     94  /// flow. It is an efficient primal-dual solution method, which
     95  /// can be viewed as the generalization of the \ref Preflow
     96  /// "preflow push-relabel" algorithm for the maximum flow problem.
    4897  ///
    49   /// \tparam Digraph The digraph type the algorithm runs on.
    50   /// \tparam LowerMap The type of the lower bound map.
    51   /// \tparam CapacityMap The type of the capacity (upper bound) map.
    52   /// \tparam CostMap The type of the cost (length) map.
    53   /// \tparam SupplyMap The type of the supply map.
     98  /// Most of the parameters of the problem (except for the digraph)
     99  /// can be given using separate functions, and the algorithm can be
     100  /// executed using the \ref run() function. If some parameters are not
     101  /// specified, then default values will be used.
    54102  ///
    55   /// \warning
    56   /// - Arc capacities and costs should be \e non-negative \e integers.
    57   /// - Supply values should be \e signed \e integers.
    58   /// - The value types of the maps should be convertible to each other.
    59   /// - \c CostMap::Value must be signed type.
     103  /// \tparam GR The digraph type the algorithm runs on.
     104  /// \tparam V The value type used for flow amounts, capacity bounds
     105  /// and supply values in the algorithm. By default it is \c int.
     106  /// \tparam C The value type used for costs and potentials in the
     107  /// algorithm. By default it is the same as \c V.
    60108  ///
    61   /// \note Arc costs are multiplied with the number of nodes during
    62   /// the algorithm so overflow problems may arise more easily than with
    63   /// other minimum cost flow algorithms.
    64   /// If it is available, <tt>long long int</tt> type is used instead of
    65   /// <tt>long int</tt> in the inside computations.
    66   ///
    67   /// \author Peter Kovacs
    68   template < typename Digraph,
    69              typename LowerMap = typename Digraph::template ArcMap<int>,
    70              typename CapacityMap = typename Digraph::template ArcMap<int>,
    71              typename CostMap = typename Digraph::template ArcMap<int>,
    72              typename SupplyMap = typename Digraph::template NodeMap<int> >
     109  /// \warning Both value types must be signed and all input data must
     110  /// be integer.
     111  /// \warning This algorithm does not support negative costs for such
     112  /// arcs that have infinite upper bound.
     113#ifdef DOXYGEN
     114  template <typename GR, typename V, typename C, typename TR>
     115#else
     116  template < typename GR, typename V = int, typename C = V,
     117             typename TR = CostScalingDefaultTraits<GR, V, C> >
     118#endif
    73119  class CostScaling
    74120  {
    75     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
    76 
    77     typedef typename CapacityMap::Value Capacity;
    78     typedef typename CostMap::Value Cost;
    79     typedef typename SupplyMap::Value Supply;
    80     typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
    81     typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
    82 
    83     typedef ResidualDigraph< const Digraph,
    84                              CapacityArcMap, CapacityArcMap > ResDigraph;
    85     typedef typename ResDigraph::Arc ResArc;
    86 
    87 #if defined __GNUC__ && !defined __STRICT_ANSI__
    88     typedef long long int LCost;
    89 #else
    90     typedef long int LCost;
    91 #endif
    92     typedef typename Digraph::template ArcMap<LCost> LargeCostMap;
    93 
    94121  public:
    95122
    96     /// The type of the flow map.
    97     typedef typename Digraph::template ArcMap<Capacity> FlowMap;
    98     /// The type of the potential map.
    99     typedef typename Digraph::template NodeMap<LCost> PotentialMap;
     123    /// The type of the digraph
     124    typedef typename TR::Digraph Digraph;
     125    /// The type of the flow amounts, capacity bounds and supply values
     126    typedef typename TR::Value Value;
     127    /// The type of the arc costs
     128    typedef typename TR::Cost Cost;
     129
     130    /// \brief The large cost type
     131    ///
     132    /// The large cost type used for internal computations.
     133    /// Using the \ref CostScalingDefaultTraits "default traits class",
     134    /// it is \c long \c long if the \c Cost type is integer,
     135    /// otherwise it is \c double.
     136    typedef typename TR::LargeCost LargeCost;
     137
     138    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
     139    typedef TR Traits;
     140
     141  public:
     142
     143    /// \brief Problem type constants for the \c run() function.
     144    ///
     145    /// Enum type containing the problem type constants that can be
     146    /// returned by the \ref run() function of the algorithm.
     147    enum ProblemType {
     148      /// The problem has no feasible solution (flow).
     149      INFEASIBLE,
     150      /// The problem has optimal solution (i.e. it is feasible and
     151      /// bounded), and the algorithm has found optimal flow and node
     152      /// potentials (primal and dual solutions).
     153      OPTIMAL,
     154      /// The digraph contains an arc of negative cost and infinite
     155      /// upper bound. It means that the objective function is unbounded
     156      /// on that arc, however note that it could actually be bounded
     157      /// over the feasible flows, but this algroithm cannot handle
     158      /// these cases.
     159      UNBOUNDED
     160    };
    100161
    101162  private:
    102163
    103     /// \brief Map adaptor class for handling residual arc costs.
    104     ///
    105     /// Map adaptor class for handling residual arc costs.
    106     template <typename Map>
    107     class ResidualCostMap : public MapBase<ResArc, typename Map::Value>
     164    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     165
     166    typedef std::vector<int> IntVector;
     167    typedef std::vector<char> BoolVector;
     168    typedef std::vector<Value> ValueVector;
     169    typedef std::vector<Cost> CostVector;
     170    typedef std::vector<LargeCost> LargeCostVector;
     171
     172  private:
     173 
     174    template <typename KT, typename VT>
     175    class VectorMap {
     176    public:
     177      typedef KT Key;
     178      typedef VT Value;
     179     
     180      VectorMap(std::vector<Value>& v) : _v(v) {}
     181     
     182      const Value& operator[](const Key& key) const {
     183        return _v[StaticDigraph::id(key)];
     184      }
     185
     186      Value& operator[](const Key& key) {
     187        return _v[StaticDigraph::id(key)];
     188      }
     189     
     190      void set(const Key& key, const Value& val) {
     191        _v[StaticDigraph::id(key)] = val;
     192      }
     193
     194    private:
     195      std::vector<Value>& _v;
     196    };
     197
     198    typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
     199    typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
     200
     201  private:
     202
     203    // Data related to the underlying digraph
     204    const GR &_graph;
     205    int _node_num;
     206    int _arc_num;
     207    int _res_node_num;
     208    int _res_arc_num;
     209    int _root;
     210
     211    // Parameters of the problem
     212    bool _have_lower;
     213    Value _sum_supply;
     214
     215    // Data structures for storing the digraph
     216    IntNodeMap _node_id;
     217    IntArcMap _arc_idf;
     218    IntArcMap _arc_idb;
     219    IntVector _first_out;
     220    BoolVector _forward;
     221    IntVector _source;
     222    IntVector _target;
     223    IntVector _reverse;
     224
     225    // Node and arc data
     226    ValueVector _lower;
     227    ValueVector _upper;
     228    CostVector _scost;
     229    ValueVector _supply;
     230
     231    ValueVector _res_cap;
     232    LargeCostVector _cost;
     233    LargeCostVector _pi;
     234    ValueVector _excess;
     235    IntVector _next_out;
     236    std::deque<int> _active_nodes;
     237
     238    // Data for scaling
     239    LargeCost _epsilon;
     240    int _alpha;
     241
     242    // Data for a StaticDigraph structure
     243    typedef std::pair<int, int> IntPair;
     244    StaticDigraph _sgr;
     245    std::vector<IntPair> _arc_vec;
     246    std::vector<LargeCost> _cost_vec;
     247    LargeCostArcMap _cost_map;
     248    LargeCostNodeMap _pi_map;
     249 
     250  public:
     251 
     252    /// \brief Constant for infinite upper bounds (capacities).
     253    ///
     254    /// Constant for infinite upper bounds (capacities).
     255    /// It is \c std::numeric_limits<Value>::infinity() if available,
     256    /// \c std::numeric_limits<Value>::max() otherwise.
     257    const Value INF;
     258
     259  public:
     260
     261    /// \name Named Template Parameters
     262    /// @{
     263
     264    template <typename T>
     265    struct SetLargeCostTraits : public Traits {
     266      typedef T LargeCost;
     267    };
     268
     269    /// \brief \ref named-templ-param "Named parameter" for setting
     270    /// \c LargeCost type.
     271    ///
     272    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
     273    /// type, which is used for internal computations in the algorithm.
     274    /// \c Cost must be convertible to \c LargeCost.
     275    template <typename T>
     276    struct SetLargeCost
     277      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
     278      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
     279    };
     280
     281    /// @}
     282
     283  public:
     284
     285    /// \brief Constructor.
     286    ///
     287    /// The constructor of the class.
     288    ///
     289    /// \param graph The digraph the algorithm runs on.
     290    CostScaling(const GR& graph) :
     291      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
     292      _cost_map(_cost_vec), _pi_map(_pi),
     293      INF(std::numeric_limits<Value>::has_infinity ?
     294          std::numeric_limits<Value>::infinity() :
     295          std::numeric_limits<Value>::max())
    108296    {
    109     private:
    110 
    111       const Map &_cost_map;
    112 
    113     public:
    114 
    115       ///\e
    116       ResidualCostMap(const Map &cost_map) :
    117         _cost_map(cost_map) {}
    118 
    119       ///\e
    120       inline typename Map::Value operator[](const ResArc &e) const {
    121         return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
    122       }
    123 
    124     }; //class ResidualCostMap
    125 
    126     /// \brief Map adaptor class for handling reduced arc costs.
    127     ///
    128     /// Map adaptor class for handling reduced arc costs.
    129     class ReducedCostMap : public MapBase<Arc, LCost>
    130     {
    131     private:
    132 
    133       const Digraph &_gr;
    134       const LargeCostMap &_cost_map;
    135       const PotentialMap &_pot_map;
    136 
    137     public:
    138 
    139       ///\e
    140       ReducedCostMap( const Digraph &gr,
    141                       const LargeCostMap &cost_map,
    142                       const PotentialMap &pot_map ) :
    143         _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
    144 
    145       ///\e
    146       inline LCost operator[](const Arc &e) const {
    147         return _cost_map[e] + _pot_map[_gr.source(e)]
    148                             - _pot_map[_gr.target(e)];
    149       }
    150 
    151     }; //class ReducedCostMap
    152 
    153   private:
    154 
    155     // The digraph the algorithm runs on
    156     const Digraph &_graph;
    157     // The original lower bound map
    158     const LowerMap *_lower;
    159     // The modified capacity map
    160     CapacityArcMap _capacity;
    161     // The original cost map
    162     const CostMap &_orig_cost;
    163     // The scaled cost map
    164     LargeCostMap _cost;
    165     // The modified supply map
    166     SupplyNodeMap _supply;
    167     bool _valid_supply;
    168 
    169     // Arc map of the current flow
    170     FlowMap *_flow;
    171     bool _local_flow;
    172     // Node map of the current potentials
    173     PotentialMap *_potential;
    174     bool _local_potential;
    175 
    176     // The residual cost map
    177     ResidualCostMap<LargeCostMap> _res_cost;
    178     // The residual digraph
    179     ResDigraph *_res_graph;
    180     // The reduced cost map
    181     ReducedCostMap *_red_cost;
    182     // The excess map
    183     SupplyNodeMap _excess;
    184     // The epsilon parameter used for cost scaling
    185     LCost _epsilon;
    186     // The scaling factor
    187     int _alpha;
    188 
    189   public:
    190 
    191     /// \brief General constructor (with lower bounds).
    192     ///
    193     /// General constructor (with lower bounds).
    194     ///
    195     /// \param digraph The digraph the algorithm runs on.
    196     /// \param lower The lower bounds of the arcs.
    197     /// \param capacity The capacities (upper bounds) of the arcs.
    198     /// \param cost The cost (length) values of the arcs.
    199     /// \param supply The supply values of the nodes (signed).
    200     CostScaling( const Digraph &digraph,
    201                  const LowerMap &lower,
    202                  const CapacityMap &capacity,
    203                  const CostMap &cost,
    204                  const SupplyMap &supply ) :
    205       _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost),
    206       _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false),
    207       _potential(NULL), _local_potential(false), _res_cost(_cost),
    208       _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
    209     {
    210       // Check the sum of supply values
    211       Supply sum = 0;
    212       for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
    213       _valid_supply = sum == 0;
     297      // Check the value types
     298      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
     299        "The flow type of CostScaling must be signed");
     300      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
     301        "The cost type of CostScaling must be signed");
     302
     303      // Resize vectors
     304      _node_num = countNodes(_graph);
     305      _arc_num = countArcs(_graph);
     306      _res_node_num = _node_num + 1;
     307      _res_arc_num = 2 * (_arc_num + _node_num);
     308      _root = _node_num;
     309
     310      _first_out.resize(_res_node_num + 1);
     311      _forward.resize(_res_arc_num);
     312      _source.resize(_res_arc_num);
     313      _target.resize(_res_arc_num);
     314      _reverse.resize(_res_arc_num);
     315
     316      _lower.resize(_res_arc_num);
     317      _upper.resize(_res_arc_num);
     318      _scost.resize(_res_arc_num);
     319      _supply.resize(_res_node_num);
    214320     
    215       for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e];
    216       for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n];
    217 
    218       // Remove non-zero lower bounds
    219       for (ArcIt e(_graph); e != INVALID; ++e) {
    220         if (lower[e] != 0) {
    221           _capacity[e] -= lower[e];
    222           _supply[_graph.source(e)] -= lower[e];
    223           _supply[_graph.target(e)] += lower[e];
    224         }
    225       }
    226     }
    227 /*
    228     /// \brief General constructor (without lower bounds).
    229     ///
    230     /// General constructor (without lower bounds).
    231     ///
    232     /// \param digraph The digraph the algorithm runs on.
    233     /// \param capacity The capacities (upper bounds) of the arcs.
    234     /// \param cost The cost (length) values of the arcs.
    235     /// \param supply The supply values of the nodes (signed).
    236     CostScaling( const Digraph &digraph,
    237                  const CapacityMap &capacity,
    238                  const CostMap &cost,
    239                  const SupplyMap &supply ) :
    240       _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
    241       _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false),
    242       _potential(NULL), _local_potential(false), _res_cost(_cost),
    243       _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
    244     {
    245       // Check the sum of supply values
    246       Supply sum = 0;
    247       for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
    248       _valid_supply = sum == 0;
    249     }
    250 
    251     /// \brief Simple constructor (with lower bounds).
    252     ///
    253     /// Simple constructor (with lower bounds).
    254     ///
    255     /// \param digraph The digraph the algorithm runs on.
    256     /// \param lower The lower bounds of the arcs.
    257     /// \param capacity The capacities (upper bounds) of the arcs.
    258     /// \param cost The cost (length) values of the arcs.
     321      _res_cap.resize(_res_arc_num);
     322      _cost.resize(_res_arc_num);
     323      _pi.resize(_res_node_num);
     324      _excess.resize(_res_node_num);
     325      _next_out.resize(_res_node_num);
     326
     327      _arc_vec.reserve(_res_arc_num);
     328      _cost_vec.reserve(_res_arc_num);
     329
     330      // Copy the graph
     331      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
     332      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     333        _node_id[n] = i;
     334      }
     335      i = 0;
     336      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     337        _first_out[i] = j;
     338        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
     339          _arc_idf[a] = j;
     340          _forward[j] = true;
     341          _source[j] = i;
     342          _target[j] = _node_id[_graph.runningNode(a)];
     343        }
     344        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
     345          _arc_idb[a] = j;
     346          _forward[j] = false;
     347          _source[j] = i;
     348          _target[j] = _node_id[_graph.runningNode(a)];
     349        }
     350        _forward[j] = false;
     351        _source[j] = i;
     352        _target[j] = _root;
     353        _reverse[j] = k;
     354        _forward[k] = true;
     355        _source[k] = _root;
     356        _target[k] = i;
     357        _reverse[k] = j;
     358        ++j; ++k;
     359      }
     360      _first_out[i] = j;
     361      _first_out[_res_node_num] = k;
     362      for (ArcIt a(_graph); a != INVALID; ++a) {
     363        int fi = _arc_idf[a];
     364        int bi = _arc_idb[a];
     365        _reverse[fi] = bi;
     366        _reverse[bi] = fi;
     367      }
     368     
     369      // Reset parameters
     370      reset();
     371    }
     372
     373    /// \name Parameters
     374    /// The parameters of the algorithm can be specified using these
     375    /// functions.
     376
     377    /// @{
     378
     379    /// \brief Set the lower bounds on the arcs.
     380    ///
     381    /// This function sets the lower bounds on the arcs.
     382    /// If it is not used before calling \ref run(), the lower bounds
     383    /// will be set to zero on all arcs.
     384    ///
     385    /// \param map An arc map storing the lower bounds.
     386    /// Its \c Value type must be convertible to the \c Value type
     387    /// of the algorithm.
     388    ///
     389    /// \return <tt>(*this)</tt>
     390    template <typename LowerMap>
     391    CostScaling& lowerMap(const LowerMap& map) {
     392      _have_lower = true;
     393      for (ArcIt a(_graph); a != INVALID; ++a) {
     394        _lower[_arc_idf[a]] = map[a];
     395        _lower[_arc_idb[a]] = map[a];
     396      }
     397      return *this;
     398    }
     399
     400    /// \brief Set the upper bounds (capacities) on the arcs.
     401    ///
     402    /// This function sets the upper bounds (capacities) on the arcs.
     403    /// If it is not used before calling \ref run(), the upper bounds
     404    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     405    /// unbounded from above on each arc).
     406    ///
     407    /// \param map An arc map storing the upper bounds.
     408    /// Its \c Value type must be convertible to the \c Value type
     409    /// of the algorithm.
     410    ///
     411    /// \return <tt>(*this)</tt>
     412    template<typename UpperMap>
     413    CostScaling& upperMap(const UpperMap& map) {
     414      for (ArcIt a(_graph); a != INVALID; ++a) {
     415        _upper[_arc_idf[a]] = map[a];
     416      }
     417      return *this;
     418    }
     419
     420    /// \brief Set the costs of the arcs.
     421    ///
     422    /// This function sets the costs of the arcs.
     423    /// If it is not used before calling \ref run(), the costs
     424    /// will be set to \c 1 on all arcs.
     425    ///
     426    /// \param map An arc map storing the costs.
     427    /// Its \c Value type must be convertible to the \c Cost type
     428    /// of the algorithm.
     429    ///
     430    /// \return <tt>(*this)</tt>
     431    template<typename CostMap>
     432    CostScaling& costMap(const CostMap& map) {
     433      for (ArcIt a(_graph); a != INVALID; ++a) {
     434        _scost[_arc_idf[a]] =  map[a];
     435        _scost[_arc_idb[a]] = -map[a];
     436      }
     437      return *this;
     438    }
     439
     440    /// \brief Set the supply values of the nodes.
     441    ///
     442    /// This function sets the supply values of the nodes.
     443    /// If neither this function nor \ref stSupply() is used before
     444    /// calling \ref run(), the supply of each node will be set to zero.
     445    ///
     446    /// \param map A node map storing the supply values.
     447    /// Its \c Value type must be convertible to the \c Value type
     448    /// of the algorithm.
     449    ///
     450    /// \return <tt>(*this)</tt>
     451    template<typename SupplyMap>
     452    CostScaling& supplyMap(const SupplyMap& map) {
     453      for (NodeIt n(_graph); n != INVALID; ++n) {
     454        _supply[_node_id[n]] = map[n];
     455      }
     456      return *this;
     457    }
     458
     459    /// \brief Set single source and target nodes and a supply value.
     460    ///
     461    /// This function sets a single source node and a single target node
     462    /// and the required flow value.
     463    /// If neither this function nor \ref supplyMap() is used before
     464    /// calling \ref run(), the supply of each node will be set to zero.
     465    ///
     466    /// Using this function has the same effect as using \ref supplyMap()
     467    /// with such a map in which \c k is assigned to \c s, \c -k is
     468    /// assigned to \c t and all other nodes have zero supply value.
     469    ///
    259470    /// \param s The source node.
    260471    /// \param t The target node.
    261     /// \param flow_value The required amount of flow from node \c s
    262     /// to node \c t (i.e. the supply of \c s and the demand of \c t).
    263     CostScaling( const Digraph &digraph,
    264                  const LowerMap &lower,
    265                  const CapacityMap &capacity,
    266                  const CostMap &cost,
    267                  Node s, Node t,
    268                  Supply flow_value ) :
    269       _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost),
    270       _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
    271       _potential(NULL), _local_potential(false), _res_cost(_cost),
    272       _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
    273     {
    274       // Remove non-zero lower bounds
    275       _supply[s] =  flow_value;
    276       _supply[t] = -flow_value;
    277       for (ArcIt e(_graph); e != INVALID; ++e) {
    278         if (lower[e] != 0) {
    279           _capacity[e] -= lower[e];
    280           _supply[_graph.source(e)] -= lower[e];
    281           _supply[_graph.target(e)] += lower[e];
    282         }
    283       }
    284       _valid_supply = true;
    285     }
    286 
    287     /// \brief Simple constructor (without lower bounds).
    288     ///
    289     /// Simple constructor (without lower bounds).
    290     ///
    291     /// \param digraph The digraph the algorithm runs on.
    292     /// \param capacity The capacities (upper bounds) of the arcs.
    293     /// \param cost The cost (length) values of the arcs.
    294     /// \param s The source node.
    295     /// \param t The target node.
    296     /// \param flow_value The required amount of flow from node \c s
    297     /// to node \c t (i.e. the supply of \c s and the demand of \c t).
    298     CostScaling( const Digraph &digraph,
    299                  const CapacityMap &capacity,
    300                  const CostMap &cost,
    301                  Node s, Node t,
    302                  Supply flow_value ) :
    303       _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
    304       _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
    305       _potential(NULL), _local_potential(false), _res_cost(_cost),
    306       _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
    307     {
    308       _supply[s] =  flow_value;
    309       _supply[t] = -flow_value;
    310       _valid_supply = true;
    311     }
    312 */
    313     /// Destructor.
    314     ~CostScaling() {
    315       if (_local_flow) delete _flow;
    316       if (_local_potential) delete _potential;
    317       delete _res_graph;
    318       delete _red_cost;
    319     }
    320 
    321     /// \brief Set the flow map.
    322     ///
    323     /// Set the flow map.
    324     ///
    325     /// \return \c (*this)
    326     CostScaling& flowMap(FlowMap &map) {
    327       if (_local_flow) {
    328         delete _flow;
    329         _local_flow = false;
    330       }
    331       _flow = &map;
     472    /// \param k The required amount of flow from node \c s to node \c t
     473    /// (i.e. the supply of \c s and the demand of \c t).
     474    ///
     475    /// \return <tt>(*this)</tt>
     476    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
     477      for (int i = 0; i != _res_node_num; ++i) {
     478        _supply[i] = 0;
     479      }
     480      _supply[_node_id[s]] =  k;
     481      _supply[_node_id[t]] = -k;
    332482      return *this;
    333483    }
    334 
    335     /// \brief Set the potential map.
    336     ///
    337     /// Set the potential map.
    338     ///
    339     /// \return \c (*this)
    340     CostScaling& potentialMap(PotentialMap &map) {
    341       if (_local_potential) {
    342         delete _potential;
    343         _local_potential = false;
    344       }
    345       _potential = &map;
    346       return *this;
    347     }
     484   
     485    /// @}
    348486
    349487    /// \name Execution control
     488    /// The algorithm can be executed using \ref run().
    350489
    351490    /// @{
     
    353492    /// \brief Run the algorithm.
    354493    ///
    355     /// Run the algorithm.
     494    /// This function runs the algorithm.
     495    /// The paramters can be specified using functions \ref lowerMap(),
     496    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     497    /// For example,
     498    /// \code
     499    ///   CostScaling<ListDigraph> cs(graph);
     500    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     501    ///     .supplyMap(sup).run();
     502    /// \endcode
     503    ///
     504    /// This function can be called more than once. All the parameters
     505    /// that have been given are kept for the next call, unless
     506    /// \ref reset() is called, thus only the modified parameters
     507    /// have to be set again. See \ref reset() for examples.
     508    /// However the underlying digraph must not be modified after this
     509    /// class have been constructed, since it copies the digraph.
    356510    ///
    357511    /// \param partial_augment By default the algorithm performs
     
    360514    /// relabel operations instead.
    361515    ///
    362     /// \return \c true if a feasible flow can be found.
    363     bool run(bool partial_augment = true) {
     516    /// \return \c INFEASIBLE if no feasible flow exists,
     517    /// \n \c OPTIMAL if the problem has optimal solution
     518    /// (i.e. it is feasible and bounded), and the algorithm has found
     519    /// optimal flow and node potentials (primal and dual solutions),
     520    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
     521    /// and infinite upper bound. It means that the objective function
     522    /// is unbounded on that arc, however note that it could actually be
     523    /// bounded over the feasible flows, but this algroithm cannot handle
     524    /// these cases.
     525    ///
     526    /// \see ProblemType
     527    ProblemType run(bool partial_augment = true) {
     528      ProblemType pt = init();
     529      if (pt != OPTIMAL) return pt;
     530      start(partial_augment);
     531      return OPTIMAL;
     532    }
     533
     534    /// \brief Reset all the parameters that have been given before.
     535    ///
     536    /// This function resets all the paramaters that have been given
     537    /// before using functions \ref lowerMap(), \ref upperMap(),
     538    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
     539    ///
     540    /// It is useful for multiple run() calls. If this function is not
     541    /// used, all the parameters given before are kept for the next
     542    /// \ref run() call.
     543    /// However the underlying digraph must not be modified after this
     544    /// class have been constructed, since it copies and extends the graph.
     545    ///
     546    /// For example,
     547    /// \code
     548    ///   CostScaling<ListDigraph> cs(graph);
     549    ///
     550    ///   // First run
     551    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     552    ///     .supplyMap(sup).run();
     553    ///
     554    ///   // Run again with modified cost map (reset() is not called,
     555    ///   // so only the cost map have to be set again)
     556    ///   cost[e] += 100;
     557    ///   cs.costMap(cost).run();
     558    ///
     559    ///   // Run again from scratch using reset()
     560    ///   // (the lower bounds will be set to zero on all arcs)
     561    ///   cs.reset();
     562    ///   cs.upperMap(capacity).costMap(cost)
     563    ///     .supplyMap(sup).run();
     564    /// \endcode
     565    ///
     566    /// \return <tt>(*this)</tt>
     567    CostScaling& reset() {
     568      for (int i = 0; i != _res_node_num; ++i) {
     569        _supply[i] = 0;
     570      }
     571      int limit = _first_out[_root];
     572      for (int j = 0; j != limit; ++j) {
     573        _lower[j] = 0;
     574        _upper[j] = INF;
     575        _scost[j] = _forward[j] ? 1 : -1;
     576      }
     577      for (int j = limit; j != _res_arc_num; ++j) {
     578        _lower[j] = 0;
     579        _upper[j] = INF;
     580        _scost[j] = 0;
     581        _scost[_reverse[j]] = 0;
     582      }     
     583      _have_lower = false;
     584      return *this;
     585    }
     586
     587    /// @}
     588
     589    /// \name Query Functions
     590    /// The results of the algorithm can be obtained using these
     591    /// functions.\n
     592    /// The \ref run() function must be called before using them.
     593
     594    /// @{
     595
     596    /// \brief Return the total cost of the found flow.
     597    ///
     598    /// This function returns the total cost of the found flow.
     599    /// Its complexity is O(e).
     600    ///
     601    /// \note The return type of the function can be specified as a
     602    /// template parameter. For example,
     603    /// \code
     604    ///   cs.totalCost<double>();
     605    /// \endcode
     606    /// It is useful if the total cost cannot be stored in the \c Cost
     607    /// type of the algorithm, which is the default return type of the
     608    /// function.
     609    ///
     610    /// \pre \ref run() must be called before using this function.
     611    template <typename Number>
     612    Number totalCost() const {
     613      Number c = 0;
     614      for (ArcIt a(_graph); a != INVALID; ++a) {
     615        int i = _arc_idb[a];
     616        c += static_cast<Number>(_res_cap[i]) *
     617             (-static_cast<Number>(_scost[i]));
     618      }
     619      return c;
     620    }
     621
     622#ifndef DOXYGEN
     623    Cost totalCost() const {
     624      return totalCost<Cost>();
     625    }
     626#endif
     627
     628    /// \brief Return the flow on the given arc.
     629    ///
     630    /// This function returns the flow on the given arc.
     631    ///
     632    /// \pre \ref run() must be called before using this function.
     633    Value flow(const Arc& a) const {
     634      return _res_cap[_arc_idb[a]];
     635    }
     636
     637    /// \brief Return the flow map (the primal solution).
     638    ///
     639    /// This function copies the flow value on each arc into the given
     640    /// map. The \c Value type of the algorithm must be convertible to
     641    /// the \c Value type of the map.
     642    ///
     643    /// \pre \ref run() must be called before using this function.
     644    template <typename FlowMap>
     645    void flowMap(FlowMap &map) const {
     646      for (ArcIt a(_graph); a != INVALID; ++a) {
     647        map.set(a, _res_cap[_arc_idb[a]]);
     648      }
     649    }
     650
     651    /// \brief Return the potential (dual value) of the given node.
     652    ///
     653    /// This function returns the potential (dual value) of the
     654    /// given node.
     655    ///
     656    /// \pre \ref run() must be called before using this function.
     657    Cost potential(const Node& n) const {
     658      return static_cast<Cost>(_pi[_node_id[n]]);
     659    }
     660
     661    /// \brief Return the potential map (the dual solution).
     662    ///
     663    /// This function copies the potential (dual value) of each node
     664    /// into the given map.
     665    /// The \c Cost type of the algorithm must be convertible to the
     666    /// \c Value type of the map.
     667    ///
     668    /// \pre \ref run() must be called before using this function.
     669    template <typename PotentialMap>
     670    void potentialMap(PotentialMap &map) const {
     671      for (NodeIt n(_graph); n != INVALID; ++n) {
     672        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
     673      }
     674    }
     675
     676    /// @}
     677
     678  private:
     679
     680    // Initialize the algorithm
     681    ProblemType init() {
     682      if (_res_node_num == 0) return INFEASIBLE;
     683
     684      // Scaling factor
     685      _alpha = 8;
     686
     687      // Check the sum of supply values
     688      _sum_supply = 0;
     689      for (int i = 0; i != _root; ++i) {
     690        _sum_supply += _supply[i];
     691      }
     692      if (_sum_supply > 0) return INFEASIBLE;
     693     
     694
     695      // Initialize vectors
     696      for (int i = 0; i != _res_node_num; ++i) {
     697        _pi[i] = 0;
     698        _excess[i] = _supply[i];
     699      }
     700     
     701      // Remove infinite upper bounds and check negative arcs
     702      const Value MAX = std::numeric_limits<Value>::max();
     703      int last_out;
     704      if (_have_lower) {
     705        for (int i = 0; i != _root; ++i) {
     706          last_out = _first_out[i+1];
     707          for (int j = _first_out[i]; j != last_out; ++j) {
     708            if (_forward[j]) {
     709              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
     710              if (c >= MAX) return UNBOUNDED;
     711              _excess[i] -= c;
     712              _excess[_target[j]] += c;
     713            }
     714          }
     715        }
     716      } else {
     717        for (int i = 0; i != _root; ++i) {
     718          last_out = _first_out[i+1];
     719          for (int j = _first_out[i]; j != last_out; ++j) {
     720            if (_forward[j] && _scost[j] < 0) {
     721              Value c = _upper[j];
     722              if (c >= MAX) return UNBOUNDED;
     723              _excess[i] -= c;
     724              _excess[_target[j]] += c;
     725            }
     726          }
     727        }
     728      }
     729      Value ex, max_cap = 0;
     730      for (int i = 0; i != _res_node_num; ++i) {
     731        ex = _excess[i];
     732        _excess[i] = 0;
     733        if (ex < 0) max_cap -= ex;
     734      }
     735      for (int j = 0; j != _res_arc_num; ++j) {
     736        if (_upper[j] >= MAX) _upper[j] = max_cap;
     737      }
     738
     739      // Initialize the large cost vector and the epsilon parameter
     740      _epsilon = 0;
     741      LargeCost lc;
     742      for (int i = 0; i != _root; ++i) {
     743        last_out = _first_out[i+1];
     744        for (int j = _first_out[i]; j != last_out; ++j) {
     745          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
     746          _cost[j] = lc;
     747          if (lc > _epsilon) _epsilon = lc;
     748        }
     749      }
     750      _epsilon /= _alpha;
     751
     752      // Initialize maps for Circulation and remove non-zero lower bounds
     753      ConstMap<Arc, Value> low(0);
     754      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
     755      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
     756      ValueArcMap cap(_graph), flow(_graph);
     757      ValueNodeMap sup(_graph);
     758      for (NodeIt n(_graph); n != INVALID; ++n) {
     759        sup[n] = _supply[_node_id[n]];
     760      }
     761      if (_have_lower) {
     762        for (ArcIt a(_graph); a != INVALID; ++a) {
     763          int j = _arc_idf[a];
     764          Value c = _lower[j];
     765          cap[a] = _upper[j] - c;
     766          sup[_graph.source(a)] -= c;
     767          sup[_graph.target(a)] += c;
     768        }
     769      } else {
     770        for (ArcIt a(_graph); a != INVALID; ++a) {
     771          cap[a] = _upper[_arc_idf[a]];
     772        }
     773      }
     774
     775      // Find a feasible flow using Circulation
     776      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
     777        circ(_graph, low, cap, sup);
     778      if (!circ.flowMap(flow).run()) return INFEASIBLE;
     779
     780      // Set residual capacities and handle GEQ supply type
     781      if (_sum_supply < 0) {
     782        for (ArcIt a(_graph); a != INVALID; ++a) {
     783          Value fa = flow[a];
     784          _res_cap[_arc_idf[a]] = cap[a] - fa;
     785          _res_cap[_arc_idb[a]] = fa;
     786          sup[_graph.source(a)] -= fa;
     787          sup[_graph.target(a)] += fa;
     788        }
     789        for (NodeIt n(_graph); n != INVALID; ++n) {
     790          _excess[_node_id[n]] = sup[n];
     791        }
     792        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
     793          int u = _target[a];
     794          int ra = _reverse[a];
     795          _res_cap[a] = -_sum_supply + 1;
     796          _res_cap[ra] = -_excess[u];
     797          _cost[a] = 0;
     798          _cost[ra] = 0;
     799          _excess[u] = 0;
     800        }
     801      } else {
     802        for (ArcIt a(_graph); a != INVALID; ++a) {
     803          Value fa = flow[a];
     804          _res_cap[_arc_idf[a]] = cap[a] - fa;
     805          _res_cap[_arc_idb[a]] = fa;
     806        }
     807        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
     808          int ra = _reverse[a];
     809          _res_cap[a] = 1;
     810          _res_cap[ra] = 0;
     811          _cost[a] = 0;
     812          _cost[ra] = 0;
     813        }
     814      }
     815     
     816      return OPTIMAL;
     817    }
     818
     819    // Execute the algorithm and transform the results
     820    void start(bool partial_augment) {
     821      // Execute the algorithm
    364822      if (partial_augment) {
    365         return init() && startPartialAugment();
     823        startPartialAugment();
    366824      } else {
    367         return init() && startPushRelabel();
    368       }
    369     }
    370 
    371     /// @}
    372 
    373     /// \name Query Functions
    374     /// The result of the algorithm can be obtained using these
    375     /// functions.\n
    376     /// \ref lemon::CostScaling::run() "run()" must be called before
    377     /// using them.
    378 
    379     /// @{
    380 
    381     /// \brief Return a const reference to the arc map storing the
    382     /// found flow.
    383     ///
    384     /// Return a const reference to the arc map storing the found flow.
    385     ///
    386     /// \pre \ref run() must be called before using this function.
    387     const FlowMap& flowMap() const {
    388       return *_flow;
    389     }
    390 
    391     /// \brief Return a const reference to the node map storing the
    392     /// found potentials (the dual solution).
    393     ///
    394     /// Return a const reference to the node map storing the found
    395     /// potentials (the dual solution).
    396     ///
    397     /// \pre \ref run() must be called before using this function.
    398     const PotentialMap& potentialMap() const {
    399       return *_potential;
    400     }
    401 
    402     /// \brief Return the flow on the given arc.
    403     ///
    404     /// Return the flow on the given arc.
    405     ///
    406     /// \pre \ref run() must be called before using this function.
    407     Capacity flow(const Arc& arc) const {
    408       return (*_flow)[arc];
    409     }
    410 
    411     /// \brief Return the potential of the given node.
    412     ///
    413     /// Return the potential of the given node.
    414     ///
    415     /// \pre \ref run() must be called before using this function.
    416     Cost potential(const Node& node) const {
    417       return (*_potential)[node];
    418     }
    419 
    420     /// \brief Return the total cost of the found flow.
    421     ///
    422     /// Return the total cost of the found flow. The complexity of the
    423     /// function is \f$ O(e) \f$.
    424     ///
    425     /// \pre \ref run() must be called before using this function.
    426     Cost totalCost() const {
    427       Cost c = 0;
    428       for (ArcIt e(_graph); e != INVALID; ++e)
    429         c += (*_flow)[e] * _orig_cost[e];
    430       return c;
    431     }
    432 
    433     /// @}
    434 
    435   private:
    436 
    437     /// Initialize the algorithm.
    438     bool init() {
    439       if (!_valid_supply) return false;
    440       // The scaling factor
    441       _alpha = 8;
    442 
    443       // Initialize flow and potential maps
    444       if (!_flow) {
    445         _flow = new FlowMap(_graph);
    446         _local_flow = true;
    447       }
    448       if (!_potential) {
    449         _potential = new PotentialMap(_graph);
    450         _local_potential = true;
    451       }
    452 
    453       _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
    454       _res_graph = new ResDigraph(_graph, _capacity, *_flow);
    455 
    456       // Initialize the scaled cost map and the epsilon parameter
    457       Cost max_cost = 0;
    458       int node_num = countNodes(_graph);
    459       for (ArcIt e(_graph); e != INVALID; ++e) {
    460         _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
    461         if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
    462       }
    463       _epsilon = max_cost * node_num;
    464 
    465       // Find a feasible flow using Circulation
    466       Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
    467                    SupplyMap >
    468         circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
    469                      _supply );
    470       return circulation.flowMap(*_flow).run();
     825        startPushRelabel();
     826      }
     827
     828      // Compute node potentials for the original costs
     829      _arc_vec.clear();
     830      _cost_vec.clear();
     831      for (int j = 0; j != _res_arc_num; ++j) {
     832        if (_res_cap[j] > 0) {
     833          _arc_vec.push_back(IntPair(_source[j], _target[j]));
     834          _cost_vec.push_back(_scost[j]);
     835        }
     836      }
     837      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
     838
     839      typename BellmanFord<StaticDigraph, LargeCostArcMap>
     840        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
     841      bf.distMap(_pi_map);
     842      bf.init(0);
     843      bf.start();
     844
     845      // Handle non-zero lower bounds
     846      if (_have_lower) {
     847        int limit = _first_out[_root];
     848        for (int j = 0; j != limit; ++j) {
     849          if (!_forward[j]) _res_cap[j] += _lower[j];
     850        }
     851      }
    471852    }
    472853
    473854    /// Execute the algorithm performing partial augmentation and
    474     /// relabel operations.
    475     bool startPartialAugment() {
     855    /// relabel operations
     856    void startPartialAugment() {
    476857      // Paramters for heuristics
    477 //      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
    478 //      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
     858      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
     859      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
    479860      // Maximum augment path length
    480861      const int MAX_PATH_LENGTH = 4;
    481862
    482       // Variables
    483       typename Digraph::template NodeMap<Arc> pred_arc(_graph);
    484       typename Digraph::template NodeMap<bool> forward(_graph);
    485       typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
    486       typename Digraph::template NodeMap<InArcIt> next_in(_graph);
    487       typename Digraph::template NodeMap<bool> next_dir(_graph);
    488       std::deque<Node> active_nodes;
    489       std::vector<Node> path_nodes;
    490 
    491 //      int node_num = countNodes(_graph);
     863      // Perform cost scaling phases
     864      IntVector pred_arc(_res_node_num);
     865      std::vector<int> path_nodes;
    492866      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    493867                                        1 : _epsilon / _alpha )
    494868      {
    495 /*
    496869        // "Early Termination" heuristic: use Bellman-Ford algorithm
    497870        // to check if the current flow is optimal
    498871        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
    499           typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
    500           ShiftCostMap shift_cost(_res_cost, 1);
    501           BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
     872          _arc_vec.clear();
     873          _cost_vec.clear();
     874          for (int j = 0; j != _res_arc_num; ++j) {
     875            if (_res_cap[j] > 0) {
     876              _arc_vec.push_back(IntPair(_source[j], _target[j]));
     877              _cost_vec.push_back(_cost[j] + 1);
     878            }
     879          }
     880          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
     881
     882          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    502883          bf.init(0);
    503884          bool done = false;
    504           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
     885          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
    505886          for (int i = 0; i < K && !done; ++i)
    506887            done = bf.processNextWeakRound();
    507888          if (done) break;
    508889        }
    509 */
     890
    510891        // Saturate arcs not satisfying the optimality condition
    511         Capacity delta;
    512         for (ArcIt e(_graph); e != INVALID; ++e) {
    513           if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
    514             delta = _capacity[e] - (*_flow)[e];
    515             _excess[_graph.source(e)] -= delta;
    516             _excess[_graph.target(e)] += delta;
    517             (*_flow)[e] = _capacity[e];
     892        for (int a = 0; a != _res_arc_num; ++a) {
     893          if (_res_cap[a] > 0 &&
     894              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
     895            Value delta = _res_cap[a];
     896            _excess[_source[a]] -= delta;
     897            _excess[_target[a]] += delta;
     898            _res_cap[a] = 0;
     899            _res_cap[_reverse[a]] += delta;
    518900          }
    519           if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
    520             _excess[_graph.target(e)] -= (*_flow)[e];
    521             _excess[_graph.source(e)] += (*_flow)[e];
    522             (*_flow)[e] = 0;
     901        }
     902       
     903        // Find active nodes (i.e. nodes with positive excess)
     904        for (int u = 0; u != _res_node_num; ++u) {
     905          if (_excess[u] > 0) _active_nodes.push_back(u);
     906        }
     907
     908        // Initialize the next arcs
     909        for (int u = 0; u != _res_node_num; ++u) {
     910          _next_out[u] = _first_out[u];
     911        }
     912
     913        // Perform partial augment and relabel operations
     914        while (true) {
     915          // Select an active node (FIFO selection)
     916          while (_active_nodes.size() > 0 &&
     917                 _excess[_active_nodes.front()] <= 0) {
     918            _active_nodes.pop_front();
    523919          }
    524         }
    525 
    526         // Find active nodes (i.e. nodes with positive excess)
    527         for (NodeIt n(_graph); n != INVALID; ++n) {
    528           if (_excess[n] > 0) active_nodes.push_back(n);
    529         }
    530 
    531         // Initialize the next arc maps
    532         for (NodeIt n(_graph); n != INVALID; ++n) {
    533           next_out[n] = OutArcIt(_graph, n);
    534           next_in[n] = InArcIt(_graph, n);
    535           next_dir[n] = true;
    536         }
    537 
    538         // Perform partial augment and relabel operations
    539         while (active_nodes.size() > 0) {
    540           // Select an active node (FIFO selection)
    541           if (_excess[active_nodes[0]] <= 0) {
    542             active_nodes.pop_front();
    543             continue;
    544           }
    545           Node start = active_nodes[0];
     920          if (_active_nodes.size() == 0) break;
     921          int start = _active_nodes.front();
    546922          path_nodes.clear();
    547923          path_nodes.push_back(start);
    548924
    549925          // Find an augmenting path from the start node
    550           Node u, tip = start;
    551           LCost min_red_cost;
    552           while ( _excess[tip] >= 0 &&
    553                   int(path_nodes.size()) <= MAX_PATH_LENGTH )
    554           {
    555             if (next_dir[tip]) {
    556               for (OutArcIt e = next_out[tip]; e != INVALID; ++e) {
    557                 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
    558                   u = _graph.target(e);
    559                   pred_arc[u] = e;
    560                   forward[u] = true;
    561                   next_out[tip] = e;
    562                   tip = u;
    563                   path_nodes.push_back(tip);
    564                   goto next_step;
    565                 }
    566               }
    567               next_dir[tip] = false;
    568             }
    569             for (InArcIt e = next_in[tip]; e != INVALID; ++e) {
    570               if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
    571                 u = _graph.source(e);
    572                 pred_arc[u] = e;
    573                 forward[u] = false;
    574                 next_in[tip] = e;
     926          int tip = start;
     927          while (_excess[tip] >= 0 &&
     928                 int(path_nodes.size()) <= MAX_PATH_LENGTH) {
     929            int u;
     930            LargeCost min_red_cost, rc;
     931            int last_out = _sum_supply < 0 ?
     932              _first_out[tip+1] : _first_out[tip+1] - 1;
     933            for (int a = _next_out[tip]; a != last_out; ++a) {
     934              if (_res_cap[a] > 0 &&
     935                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
     936                u = _target[a];
     937                pred_arc[u] = a;
     938                _next_out[tip] = a;
    575939                tip = u;
    576940                path_nodes.push_back(tip);
     
    580944
    581945            // Relabel tip node
    582             min_red_cost = std::numeric_limits<LCost>::max() / 2;
    583             for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) {
    584               if ( _capacity[oe] - (*_flow)[oe] > 0 &&
    585                    (*_red_cost)[oe] < min_red_cost )
    586                 min_red_cost = (*_red_cost)[oe];
     946            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
     947            for (int a = _first_out[tip]; a != last_out; ++a) {
     948              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
     949              if (_res_cap[a] > 0 && rc < min_red_cost) {
     950                min_red_cost = rc;
     951              }
    587952            }
    588             for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) {
    589               if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
    590                 min_red_cost = -(*_red_cost)[ie];
    591             }
    592             (*_potential)[tip] -= min_red_cost + _epsilon;
    593 
    594             // Reset the next arc maps
    595             next_out[tip] = OutArcIt(_graph, tip);
    596             next_in[tip] = InArcIt(_graph, tip);
    597             next_dir[tip] = true;
     953            _pi[tip] -= min_red_cost + _epsilon;
     954
     955            // Reset the next arc of tip
     956            _next_out[tip] = _first_out[tip];
    598957
    599958            // Step back
    600959            if (tip != start) {
    601960              path_nodes.pop_back();
    602               tip = path_nodes[path_nodes.size()-1];
     961              tip = path_nodes.back();
    603962            }
    604963
    605           next_step:
    606             continue;
     964          next_step: ;
    607965          }
    608966
    609967          // Augment along the found path (as much flow as possible)
    610           Capacity delta;
     968          Value delta;
     969          int u, v = path_nodes.front(), pa;
    611970          for (int i = 1; i < int(path_nodes.size()); ++i) {
    612             u = path_nodes[i];
    613             delta = forward[u] ?
    614               _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] :
    615               (*_flow)[pred_arc[u]];
    616             delta = std::min(delta, _excess[path_nodes[i-1]]);
    617             (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta;
    618             _excess[path_nodes[i-1]] -= delta;
    619             _excess[u] += delta;
    620             if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
     971            u = v;
     972            v = path_nodes[i];
     973            pa = pred_arc[v];
     974            delta = std::min(_res_cap[pa], _excess[u]);
     975            _res_cap[pa] -= delta;
     976            _res_cap[_reverse[pa]] += delta;
     977            _excess[u] -= delta;
     978            _excess[v] += delta;
     979            if (_excess[v] > 0 && _excess[v] <= delta)
     980              _active_nodes.push_back(v);
    621981          }
    622982        }
    623983      }
    624 
    625       // Compute node potentials for the original costs
    626       ResidualCostMap<CostMap> res_cost(_orig_cost);
    627       BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
    628         bf(*_res_graph, res_cost);
    629       bf.init(0); bf.start();
    630       for (NodeIt n(_graph); n != INVALID; ++n)
    631         (*_potential)[n] = bf.dist(n);
    632 
    633       // Handle non-zero lower bounds
    634       if (_lower) {
    635         for (ArcIt e(_graph); e != INVALID; ++e)
    636           (*_flow)[e] += (*_lower)[e];
    637       }
    638       return true;
    639     }
    640 
    641     /// Execute the algorithm performing push and relabel operations.
    642     bool startPushRelabel() {
     984    }
     985
     986    /// Execute the algorithm performing push and relabel operations
     987    void startPushRelabel() {
    643988      // Paramters for heuristics
    644 //      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
    645 //      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
    646 
    647       typename Digraph::template NodeMap<bool> hyper(_graph, false);
    648       typename Digraph::template NodeMap<Arc> pred_arc(_graph);
    649       typename Digraph::template NodeMap<bool> forward(_graph);
    650       typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
    651       typename Digraph::template NodeMap<InArcIt> next_in(_graph);
    652       typename Digraph::template NodeMap<bool> next_dir(_graph);
    653       std::deque<Node> active_nodes;
    654 
    655 //      int node_num = countNodes(_graph);
     989      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
     990      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
     991
     992      // Perform cost scaling phases
     993      BoolVector hyper(_res_node_num, false);
    656994      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    657995                                        1 : _epsilon / _alpha )
    658996      {
    659 /*
    660997        // "Early Termination" heuristic: use Bellman-Ford algorithm
    661998        // to check if the current flow is optimal
    662999        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
    663           typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
    664           ShiftCostMap shift_cost(_res_cost, 1);
    665           BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
     1000          _arc_vec.clear();
     1001          _cost_vec.clear();
     1002          for (int j = 0; j != _res_arc_num; ++j) {
     1003            if (_res_cap[j] > 0) {
     1004              _arc_vec.push_back(IntPair(_source[j], _target[j]));
     1005              _cost_vec.push_back(_cost[j] + 1);
     1006            }
     1007          }
     1008          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
     1009
     1010          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    6661011          bf.init(0);
    6671012          bool done = false;
    668           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
     1013          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
    6691014          for (int i = 0; i < K && !done; ++i)
    6701015            done = bf.processNextWeakRound();
    6711016          if (done) break;
    6721017        }
    673 */
    6741018
    6751019        // Saturate arcs not satisfying the optimality condition
    676         Capacity delta;
    677         for (ArcIt e(_graph); e != INVALID; ++e) {
    678           if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
    679             delta = _capacity[e] - (*_flow)[e];
    680             _excess[_graph.source(e)] -= delta;
    681             _excess[_graph.target(e)] += delta;
    682             (*_flow)[e] = _capacity[e];
     1020        for (int a = 0; a != _res_arc_num; ++a) {
     1021          if (_res_cap[a] > 0 &&
     1022              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
     1023            Value delta = _res_cap[a];
     1024            _excess[_source[a]] -= delta;
     1025            _excess[_target[a]] += delta;
     1026            _res_cap[a] = 0;
     1027            _res_cap[_reverse[a]] += delta;
    6831028          }
    684           if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
    685             _excess[_graph.target(e)] -= (*_flow)[e];
    686             _excess[_graph.source(e)] += (*_flow)[e];
    687             (*_flow)[e] = 0;
    688           }
    6891029        }
    6901030
    6911031        // Find active nodes (i.e. nodes with positive excess)
    692         for (NodeIt n(_graph); n != INVALID; ++n) {
    693           if (_excess[n] > 0) active_nodes.push_back(n);
    694         }
    695 
    696         // Initialize the next arc maps
    697         for (NodeIt n(_graph); n != INVALID; ++n) {
    698           next_out[n] = OutArcIt(_graph, n);
    699           next_in[n] = InArcIt(_graph, n);
    700           next_dir[n] = true;
     1032        for (int u = 0; u != _res_node_num; ++u) {
     1033          if (_excess[u] > 0) _active_nodes.push_back(u);
     1034        }
     1035
     1036        // Initialize the next arcs
     1037        for (int u = 0; u != _res_node_num; ++u) {
     1038          _next_out[u] = _first_out[u];
    7011039        }
    7021040
    7031041        // Perform push and relabel operations
    704         while (active_nodes.size() > 0) {
     1042        while (_active_nodes.size() > 0) {
     1043          LargeCost min_red_cost, rc;
     1044          Value delta;
     1045          int n, t, a, last_out = _res_arc_num;
     1046
    7051047          // Select an active node (FIFO selection)
    706           Node n = active_nodes[0], t;
    707           bool relabel_enabled = true;
     1048        next_node:
     1049          n = _active_nodes.front();
     1050          last_out = _sum_supply < 0 ?
     1051            _first_out[n+1] : _first_out[n+1] - 1;
    7081052
    7091053          // Perform push operations if there are admissible arcs
    710           if (_excess[n] > 0 && next_dir[n]) {
    711             OutArcIt e = next_out[n];
    712             for ( ; e != INVALID; ++e) {
    713               if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
    714                 delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
    715                 t = _graph.target(e);
     1054          if (_excess[n] > 0) {
     1055            for (a = _next_out[n]; a != last_out; ++a) {
     1056              if (_res_cap[a] > 0 &&
     1057                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
     1058                delta = std::min(_res_cap[a], _excess[n]);
     1059                t = _target[a];
    7161060
    7171061                // Push-look-ahead heuristic
    718                 Capacity ahead = -_excess[t];
    719                 for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
    720                   if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
    721                     ahead += _capacity[oe] - (*_flow)[oe];
    722                 }
    723                 for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
    724                   if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
    725                     ahead += (*_flow)[ie];
     1062                Value ahead = -_excess[t];
     1063                int last_out_t = _sum_supply < 0 ?
     1064                  _first_out[t+1] : _first_out[t+1] - 1;
     1065                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
     1066                  if (_res_cap[ta] > 0 &&
     1067                      _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
     1068                    ahead += _res_cap[ta];
     1069                  if (ahead >= delta) break;
    7261070                }
    7271071                if (ahead < 0) ahead = 0;
     
    7291073                // Push flow along the arc
    7301074                if (ahead < delta) {
    731                   (*_flow)[e] += ahead;
     1075                  _res_cap[a] -= ahead;
     1076                  _res_cap[_reverse[a]] += ahead;
    7321077                  _excess[n] -= ahead;
    7331078                  _excess[t] += ahead;
    734                   active_nodes.push_front(t);
     1079                  _active_nodes.push_front(t);
    7351080                  hyper[t] = true;
    736                   relabel_enabled = false;
    737                   break;
     1081                  _next_out[n] = a;
     1082                  goto next_node;
    7381083                } else {
    739                   (*_flow)[e] += delta;
     1084                  _res_cap[a] -= delta;
     1085                  _res_cap[_reverse[a]] += delta;
    7401086                  _excess[n] -= delta;
    7411087                  _excess[t] += delta;
    7421088                  if (_excess[t] > 0 && _excess[t] <= delta)
    743                     active_nodes.push_back(t);
     1089                    _active_nodes.push_back(t);
    7441090                }
    7451091
    746                 if (_excess[n] == 0) break;
     1092                if (_excess[n] == 0) {
     1093                  _next_out[n] = a;
     1094                  goto remove_nodes;
     1095                }
    7471096              }
    7481097            }
    749             if (e != INVALID) {
    750               next_out[n] = e;
    751             } else {
    752               next_dir[n] = false;
    753             }
     1098            _next_out[n] = a;
    7541099          }
    7551100
    756           if (_excess[n] > 0 && !next_dir[n]) {
    757             InArcIt e = next_in[n];
    758             for ( ; e != INVALID; ++e) {
    759               if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
    760                 delta = std::min((*_flow)[e], _excess[n]);
    761                 t = _graph.source(e);
    762 
    763                 // Push-look-ahead heuristic
    764                 Capacity ahead = -_excess[t];
    765                 for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
    766                   if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
    767                     ahead += _capacity[oe] - (*_flow)[oe];
    768                 }
    769                 for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
    770                   if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
    771                     ahead += (*_flow)[ie];
    772                 }
    773                 if (ahead < 0) ahead = 0;
    774 
    775                 // Push flow along the arc
    776                 if (ahead < delta) {
    777                   (*_flow)[e] -= ahead;
    778                   _excess[n] -= ahead;
    779                   _excess[t] += ahead;
    780                   active_nodes.push_front(t);
    781                   hyper[t] = true;
    782                   relabel_enabled = false;
    783                   break;
    784                 } else {
    785                   (*_flow)[e] -= delta;
    786                   _excess[n] -= delta;
    787                   _excess[t] += delta;
    788                   if (_excess[t] > 0 && _excess[t] <= delta)
    789                     active_nodes.push_back(t);
    790                 }
    791 
    792                 if (_excess[n] == 0) break;
     1101          // Relabel the node if it is still active (or hyper)
     1102          if (_excess[n] > 0 || hyper[n]) {
     1103            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
     1104            for (int a = _first_out[n]; a != last_out; ++a) {
     1105              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
     1106              if (_res_cap[a] > 0 && rc < min_red_cost) {
     1107                min_red_cost = rc;
    7931108              }
    7941109            }
    795             next_in[n] = e;
     1110            _pi[n] -= min_red_cost + _epsilon;
     1111            hyper[n] = false;
     1112
     1113            // Reset the next arc
     1114            _next_out[n] = _first_out[n];
    7961115          }
    797 
    798           // Relabel the node if it is still active (or hyper)
    799           if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
    800             LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
    801             for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) {
    802               if ( _capacity[oe] - (*_flow)[oe] > 0 &&
    803                    (*_red_cost)[oe] < min_red_cost )
    804                 min_red_cost = (*_red_cost)[oe];
    805             }
    806             for (InArcIt ie(_graph, n); ie != INVALID; ++ie) {
    807               if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
    808                 min_red_cost = -(*_red_cost)[ie];
    809             }
    810             (*_potential)[n] -= min_red_cost + _epsilon;
    811             hyper[n] = false;
    812 
    813             // Reset the next arc maps
    814             next_out[n] = OutArcIt(_graph, n);
    815             next_in[n] = InArcIt(_graph, n);
    816             next_dir[n] = true;
     1116       
     1117          // Remove nodes that are not active nor hyper
     1118        remove_nodes:
     1119          while ( _active_nodes.size() > 0 &&
     1120                  _excess[_active_nodes.front()] <= 0 &&
     1121                  !hyper[_active_nodes.front()] ) {
     1122            _active_nodes.pop_front();
    8171123          }
    818 
    819           // Remove nodes that are not active nor hyper
    820           while ( active_nodes.size() > 0 &&
    821                   _excess[active_nodes[0]] <= 0 &&
    822                   !hyper[active_nodes[0]] ) {
    823             active_nodes.pop_front();
    824           }
    825         }
    826       }
    827 
    828       // Compute node potentials for the original costs
    829       ResidualCostMap<CostMap> res_cost(_orig_cost);
    830       BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
    831         bf(*_res_graph, res_cost);
    832       bf.init(0); bf.start();
    833       for (NodeIt n(_graph); n != INVALID; ++n)
    834         (*_potential)[n] = bf.dist(n);
    835 
    836       // Handle non-zero lower bounds
    837       if (_lower) {
    838         for (ArcIt e(_graph); e != INVALID; ++e)
    839           (*_flow)[e] += (*_lower)[e];
    840       }
    841       return true;
     1124        }
     1125      }
    8421126    }
    8431127
Note: See TracChangeset for help on using the changeset viewer.