lemon/capacity_scaling.h
author alpar
Tue, 05 Jun 2007 17:30:52 +0000
changeset 2452 24887f3a35ec
child 2457 8c791ee69a45
permissions -rw-r--r--
Update the demo file of Circulation
deba@2440
     1
/* -*- C++ -*-
deba@2440
     2
 *
deba@2440
     3
 * This file is a part of LEMON, a generic C++ optimization library
deba@2440
     4
 *
deba@2440
     5
 * Copyright (C) 2003-2007
deba@2440
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@2440
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@2440
     8
 *
deba@2440
     9
 * Permission to use, modify and distribute this software is granted
deba@2440
    10
 * provided that this copyright notice appears in all copies. For
deba@2440
    11
 * precise terms see the accompanying LICENSE file.
deba@2440
    12
 *
deba@2440
    13
 * This software is provided "AS IS" with no warranty of any kind,
deba@2440
    14
 * express or implied, and with no claim as to its suitability for any
deba@2440
    15
 * purpose.
deba@2440
    16
 *
deba@2440
    17
 */
deba@2440
    18
deba@2440
    19
#ifndef LEMON_CAPACITY_SCALING_H
deba@2440
    20
#define LEMON_CAPACITY_SCALING_H
deba@2440
    21
deba@2440
    22
/// \ingroup min_cost_flow
deba@2440
    23
///
deba@2440
    24
/// \file
deba@2440
    25
/// \brief The capacity scaling algorithm for finding a minimum cost
deba@2440
    26
/// flow.
deba@2440
    27
deba@2440
    28
#include <vector>
deba@2440
    29
#include <lemon/dijkstra.h>
deba@2440
    30
#include <lemon/graph_adaptor.h>
deba@2440
    31
deba@2440
    32
#define WITH_SCALING
deba@2440
    33
deba@2440
    34
namespace lemon {
deba@2440
    35
deba@2440
    36
  /// \addtogroup min_cost_flow
deba@2440
    37
  /// @{
deba@2440
    38
deba@2440
    39
  /// \brief Implementation of the capacity scaling version of the
deba@2440
    40
  /// succesive shortest path algorithm for finding a minimum cost flow.
deba@2440
    41
  ///
deba@2440
    42
  /// \ref lemon::CapacityScaling "CapacityScaling" implements a
deba@2440
    43
  /// capacity scaling algorithm for finding a minimum cost flow.
deba@2440
    44
  ///
deba@2440
    45
  /// \param Graph The directed graph type the algorithm runs on.
deba@2440
    46
  /// \param LowerMap The type of the lower bound map.
deba@2440
    47
  /// \param CapacityMap The type of the capacity (upper bound) map.
deba@2440
    48
  /// \param CostMap The type of the cost (length) map.
deba@2440
    49
  /// \param SupplyMap The type of the supply map.
deba@2440
    50
  ///
deba@2440
    51
  /// \warning
deba@2440
    52
  /// - Edge capacities and costs should be nonnegative integers.
deba@2440
    53
  ///	However \c CostMap::Value should be signed type.
deba@2440
    54
  /// - Supply values should be integers.
deba@2440
    55
  /// - \c LowerMap::Value must be convertible to
deba@2440
    56
  ///	\c CapacityMap::Value and \c CapacityMap::Value must be
deba@2440
    57
  ///	convertible to \c SupplyMap::Value.
deba@2440
    58
  ///
deba@2440
    59
  /// \author Peter Kovacs
deba@2440
    60
deba@2440
    61
template < typename Graph,
deba@2440
    62
	   typename LowerMap = typename Graph::template EdgeMap<int>,
deba@2440
    63
	   typename CapacityMap = LowerMap,
deba@2440
    64
	   typename CostMap = typename Graph::template EdgeMap<int>,
deba@2440
    65
	   typename SupplyMap = typename Graph::template NodeMap
deba@2440
    66
				<typename CapacityMap::Value> >
deba@2440
    67
  class CapacityScaling
deba@2440
    68
  {
deba@2440
    69
    typedef typename Graph::Node Node;
deba@2440
    70
    typedef typename Graph::NodeIt NodeIt;
deba@2440
    71
    typedef typename Graph::Edge Edge;
deba@2440
    72
    typedef typename Graph::EdgeIt EdgeIt;
deba@2440
    73
    typedef typename Graph::InEdgeIt InEdgeIt;
deba@2440
    74
    typedef typename Graph::OutEdgeIt OutEdgeIt;
deba@2440
    75
deba@2440
    76
    typedef typename LowerMap::Value Lower;
deba@2440
    77
    typedef typename CapacityMap::Value Capacity;
deba@2440
    78
    typedef typename CostMap::Value Cost;
deba@2440
    79
    typedef typename SupplyMap::Value Supply;
deba@2440
    80
    typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
deba@2440
    81
    typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
deba@2440
    82
deba@2440
    83
    typedef ResGraphAdaptor< const Graph, Capacity,
deba@2440
    84
			     CapacityRefMap, CapacityRefMap > ResGraph;
deba@2440
    85
    typedef typename ResGraph::Node ResNode;
deba@2440
    86
    typedef typename ResGraph::NodeIt ResNodeIt;
deba@2440
    87
    typedef typename ResGraph::Edge ResEdge;
deba@2440
    88
    typedef typename ResGraph::EdgeIt ResEdgeIt;
deba@2440
    89
deba@2440
    90
  public:
deba@2440
    91
deba@2440
    92
    /// \brief The type of the flow map.
deba@2440
    93
    typedef CapacityRefMap FlowMap;
deba@2440
    94
    /// \brief The type of the potential map.
deba@2440
    95
    typedef typename Graph::template NodeMap<Cost> PotentialMap;
deba@2440
    96
deba@2440
    97
  protected:
deba@2440
    98
deba@2440
    99
    /// \brief Map adaptor class for handling reduced edge costs.
deba@2440
   100
    class ReducedCostMap : public MapBase<ResEdge, Cost>
deba@2440
   101
    {
deba@2440
   102
    private:
deba@2440
   103
deba@2440
   104
      const ResGraph &gr;
deba@2440
   105
      const CostMap &cost_map;
deba@2440
   106
      const PotentialMap &pot_map;
deba@2440
   107
deba@2440
   108
    public:
deba@2440
   109
deba@2440
   110
      typedef typename MapBase<ResEdge, Cost>::Value Value;
deba@2440
   111
      typedef typename MapBase<ResEdge, Cost>::Key Key;
deba@2440
   112
deba@2440
   113
      ReducedCostMap( const ResGraph &_gr,
deba@2440
   114
		      const CostMap &_cost,
deba@2440
   115
		      const PotentialMap &_pot ) :
deba@2440
   116
	gr(_gr), cost_map(_cost), pot_map(_pot) {}
deba@2440
   117
deba@2440
   118
      Value operator[](const Key &e) const {
deba@2440
   119
	return ResGraph::forward(e) ?
deba@2440
   120
	   cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :
deba@2440
   121
	  -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
deba@2440
   122
      }
deba@2440
   123
deba@2440
   124
    }; //class ReducedCostMap
deba@2440
   125
deba@2440
   126
    /// \brief Map adaptor for \ref lemon::Dijkstra "Dijkstra" class to
deba@2440
   127
    /// update node potentials.
deba@2440
   128
    class PotentialUpdateMap : public MapBase<ResNode, Cost>
deba@2440
   129
    {
deba@2440
   130
    private:
deba@2440
   131
deba@2440
   132
      PotentialMap *pot;
deba@2440
   133
      typedef std::pair<ResNode, Cost> Pair;
deba@2440
   134
      std::vector<Pair> data;
deba@2440
   135
deba@2440
   136
    public:
deba@2440
   137
deba@2440
   138
      typedef typename MapBase<ResNode, Cost>::Value Value;
deba@2440
   139
      typedef typename MapBase<ResNode, Cost>::Key Key;
deba@2440
   140
deba@2440
   141
      void potentialMap(PotentialMap &_pot) {
deba@2440
   142
	pot = &_pot;
deba@2440
   143
      }
deba@2440
   144
deba@2440
   145
      void init() {
deba@2440
   146
	data.clear();
deba@2440
   147
      }
deba@2440
   148
deba@2440
   149
      void set(const Key &n, const Value &v) {
deba@2440
   150
	data.push_back(Pair(n, v));
deba@2440
   151
      }
deba@2440
   152
deba@2440
   153
      void update() {
deba@2440
   154
	Cost val = data[data.size()-1].second;
deba@2440
   155
	for (int i = 0; i < data.size()-1; ++i)
deba@2440
   156
	  (*pot)[data[i].first] += val - data[i].second;
deba@2440
   157
      }
deba@2440
   158
deba@2440
   159
    }; //class PotentialUpdateMap
deba@2440
   160
deba@2440
   161
#ifdef WITH_SCALING
deba@2440
   162
    /// \brief Map adaptor class for identifing deficit nodes.
deba@2440
   163
    class DeficitBoolMap : public MapBase<ResNode, bool>
deba@2440
   164
    {
deba@2440
   165
    private:
deba@2440
   166
deba@2440
   167
      const SupplyRefMap &imb;
deba@2440
   168
      const Capacity &delta;
deba@2440
   169
deba@2440
   170
    public:
deba@2440
   171
deba@2440
   172
      DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :
deba@2440
   173
	imb(_imb), delta(_delta) {}
deba@2440
   174
deba@2440
   175
      bool operator[](const ResNode &n) const {
deba@2440
   176
	return imb[n] <= -delta;
deba@2440
   177
      }
deba@2440
   178
deba@2440
   179
    }; //class DeficitBoolMap
deba@2440
   180
deba@2440
   181
    /// \brief Map adaptor class for filtering edges with at least
deba@2440
   182
    /// \c delta residual capacity
deba@2440
   183
    class DeltaFilterMap : public MapBase<ResEdge, bool>
deba@2440
   184
    {
deba@2440
   185
    private:
deba@2440
   186
deba@2440
   187
      const ResGraph &gr;
deba@2440
   188
      const Capacity &delta;
deba@2440
   189
deba@2440
   190
    public:
deba@2440
   191
deba@2440
   192
      typedef typename MapBase<ResEdge, Cost>::Value Value;
deba@2440
   193
      typedef typename MapBase<ResEdge, Cost>::Key Key;
deba@2440
   194
deba@2440
   195
      DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :
deba@2440
   196
	gr(_gr), delta(_delta) {}
deba@2440
   197
deba@2440
   198
      Value operator[](const Key &e) const {
deba@2440
   199
	return gr.rescap(e) >= delta;
deba@2440
   200
      }
deba@2440
   201
deba@2440
   202
    }; //class DeltaFilterMap
deba@2440
   203
deba@2440
   204
    typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>
deba@2440
   205
      DeltaResGraph;
deba@2440
   206
deba@2440
   207
    /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
deba@2440
   208
    class ResDijkstraTraits :
deba@2440
   209
      public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>
deba@2440
   210
    {
deba@2440
   211
    public:
deba@2440
   212
deba@2440
   213
      typedef PotentialUpdateMap DistMap;
deba@2440
   214
deba@2440
   215
      static DistMap *createDistMap(const DeltaResGraph&) {
deba@2440
   216
	return new DistMap();
deba@2440
   217
      }
deba@2440
   218
deba@2440
   219
    }; //class ResDijkstraTraits
deba@2440
   220
deba@2440
   221
#else //WITHOUT_CAPACITY_SCALING
deba@2440
   222
    /// \brief Map adaptor class for identifing deficit nodes.
deba@2440
   223
    class DeficitBoolMap : public MapBase<ResNode, bool>
deba@2440
   224
    {
deba@2440
   225
    private:
deba@2440
   226
deba@2440
   227
      const SupplyRefMap &imb;
deba@2440
   228
deba@2440
   229
    public:
deba@2440
   230
deba@2440
   231
      DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}
deba@2440
   232
deba@2440
   233
      bool operator[](const ResNode &n) const {
deba@2440
   234
	return imb[n] < 0;
deba@2440
   235
      }
deba@2440
   236
deba@2440
   237
    }; //class DeficitBoolMap
deba@2440
   238
deba@2440
   239
    /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
deba@2440
   240
    class ResDijkstraTraits :
deba@2440
   241
      public DijkstraDefaultTraits<ResGraph, ReducedCostMap>
deba@2440
   242
    {
deba@2440
   243
    public:
deba@2440
   244
deba@2440
   245
      typedef PotentialUpdateMap DistMap;
deba@2440
   246
deba@2440
   247
      static DistMap *createDistMap(const ResGraph&) {
deba@2440
   248
	return new DistMap();
deba@2440
   249
      }
deba@2440
   250
deba@2440
   251
    }; //class ResDijkstraTraits
deba@2440
   252
#endif
deba@2440
   253
deba@2440
   254
  protected:
deba@2440
   255
deba@2440
   256
    /// \brief The directed graph the algorithm runs on.
deba@2440
   257
    const Graph &graph;
deba@2440
   258
    /// \brief The original lower bound map.
deba@2440
   259
    const LowerMap *lower;
deba@2440
   260
    /// \brief The modified capacity map.
deba@2440
   261
    CapacityRefMap capacity;
deba@2440
   262
    /// \brief The cost map.
deba@2440
   263
    const CostMap &cost;
deba@2440
   264
    /// \brief The modified supply map.
deba@2440
   265
    SupplyRefMap supply;
deba@2440
   266
    /// \brief The sum of supply values equals zero.
deba@2440
   267
    bool valid_supply;
deba@2440
   268
deba@2440
   269
    /// \brief The edge map of the current flow.
deba@2440
   270
    FlowMap flow;
deba@2440
   271
    /// \brief The potential node map.
deba@2440
   272
    PotentialMap potential;
deba@2440
   273
    /// \brief The residual graph.
deba@2440
   274
    ResGraph res_graph;
deba@2440
   275
    /// \brief The reduced cost map.
deba@2440
   276
    ReducedCostMap red_cost;
deba@2440
   277
deba@2440
   278
    /// \brief The imbalance map.
deba@2440
   279
    SupplyRefMap imbalance;
deba@2440
   280
    /// \brief The excess nodes.
deba@2440
   281
    std::vector<ResNode> excess_nodes;
deba@2440
   282
    /// \brief The index of the next excess node.
deba@2440
   283
    int next_node;
deba@2440
   284
deba@2440
   285
#ifdef WITH_SCALING
deba@2440
   286
    typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>
deba@2440
   287
      ResDijkstra;
deba@2440
   288
    /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
deba@2440
   289
    /// shortest paths in the residual graph with respect to the
deba@2440
   290
    /// reduced edge costs.
deba@2440
   291
    ResDijkstra dijkstra;
deba@2440
   292
deba@2440
   293
    /// \brief The delta parameter used for capacity scaling.
deba@2440
   294
    Capacity delta;
deba@2440
   295
    /// \brief Edge filter map.
deba@2440
   296
    DeltaFilterMap delta_filter;
deba@2440
   297
    /// \brief The delta residual graph.
deba@2440
   298
    DeltaResGraph dres_graph;
deba@2440
   299
    /// \brief Map for identifing deficit nodes.
deba@2440
   300
    DeficitBoolMap delta_deficit;
deba@2440
   301
deba@2440
   302
#else //WITHOUT_CAPACITY_SCALING
deba@2440
   303
    typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
deba@2440
   304
      ResDijkstra;
deba@2440
   305
    /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
deba@2440
   306
    /// shortest paths in the residual graph with respect to the
deba@2440
   307
    /// reduced edge costs.
deba@2440
   308
    ResDijkstra dijkstra;
deba@2440
   309
    /// \brief Map for identifing deficit nodes.
deba@2440
   310
    DeficitBoolMap has_deficit;
deba@2440
   311
#endif
deba@2440
   312
deba@2440
   313
    /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.
deba@2440
   314
    typename ResDijkstra::PredMap pred;
deba@2440
   315
    /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to
deba@2440
   316
    /// update node potentials.
deba@2440
   317
    PotentialUpdateMap updater;
deba@2440
   318
deba@2440
   319
  public :
deba@2440
   320
deba@2440
   321
    /// \brief General constructor of the class (with lower bounds).
deba@2440
   322
    ///
deba@2440
   323
    /// General constructor of the class (with lower bounds).
deba@2440
   324
    ///
deba@2440
   325
    /// \param _graph The directed graph the algorithm runs on.
deba@2440
   326
    /// \param _lower The lower bounds of the edges.
deba@2440
   327
    /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440
   328
    /// \param _cost The cost (length) values of the edges.
deba@2440
   329
    /// \param _supply The supply values of the nodes (signed).
deba@2440
   330
    CapacityScaling( const Graph &_graph,
deba@2440
   331
		     const LowerMap &_lower,
deba@2440
   332
		     const CapacityMap &_capacity,
deba@2440
   333
		     const CostMap &_cost,
deba@2440
   334
		     const SupplyMap &_supply ) :
deba@2440
   335
      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
deba@2440
   336
      supply(_graph), flow(_graph, 0), potential(_graph, 0),
deba@2440
   337
      res_graph(_graph, capacity, flow),
deba@2440
   338
      red_cost(res_graph, cost, potential), imbalance(_graph),
deba@2440
   339
#ifdef WITH_SCALING
deba@2440
   340
      delta(0), delta_filter(res_graph, delta),
deba@2440
   341
      dres_graph(res_graph, delta_filter),
deba@2440
   342
      dijkstra(dres_graph, red_cost), pred(dres_graph),
deba@2440
   343
      delta_deficit(imbalance, delta)
deba@2440
   344
#else //WITHOUT_CAPACITY_SCALING
deba@2440
   345
      dijkstra(res_graph, red_cost), pred(res_graph),
deba@2440
   346
      has_deficit(imbalance)
deba@2440
   347
#endif
deba@2440
   348
    {
deba@2440
   349
      // Removing nonzero lower bounds
deba@2440
   350
      capacity = subMap(_capacity, _lower);
deba@2440
   351
      Supply sum = 0;
deba@2440
   352
      for (NodeIt n(graph); n != INVALID; ++n) {
deba@2440
   353
	Supply s = _supply[n];
deba@2440
   354
	for (InEdgeIt e(graph, n); e != INVALID; ++e)
deba@2440
   355
	  s += _lower[e];
deba@2440
   356
	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
deba@2440
   357
	  s -= _lower[e];
deba@2440
   358
	supply[n] = imbalance[n] = s;
deba@2440
   359
	sum += s;
deba@2440
   360
      }
deba@2440
   361
      valid_supply = sum == 0;
deba@2440
   362
    }
deba@2440
   363
deba@2440
   364
    /// \brief General constructor of the class (without lower bounds).
deba@2440
   365
    ///
deba@2440
   366
    /// General constructor of the class (without lower bounds).
deba@2440
   367
    ///
deba@2440
   368
    /// \param _graph The directed graph the algorithm runs on.
deba@2440
   369
    /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440
   370
    /// \param _cost The cost (length) values of the edges.
deba@2440
   371
    /// \param _supply The supply values of the nodes (signed).
deba@2440
   372
    CapacityScaling( const Graph &_graph,
deba@2440
   373
		     const CapacityMap &_capacity,
deba@2440
   374
		     const CostMap &_cost,
deba@2440
   375
		     const SupplyMap &_supply ) :
deba@2440
   376
      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
deba@2440
   377
      supply(_supply), flow(_graph, 0), potential(_graph, 0),
deba@2440
   378
      res_graph(_graph, capacity, flow),
deba@2440
   379
      red_cost(res_graph, cost, potential), imbalance(_graph),
deba@2440
   380
#ifdef WITH_SCALING
deba@2440
   381
      delta(0), delta_filter(res_graph, delta),
deba@2440
   382
      dres_graph(res_graph, delta_filter),
deba@2440
   383
      dijkstra(dres_graph, red_cost), pred(dres_graph),
deba@2440
   384
      delta_deficit(imbalance, delta)
deba@2440
   385
#else //WITHOUT_CAPACITY_SCALING
deba@2440
   386
      dijkstra(res_graph, red_cost), pred(res_graph),
deba@2440
   387
      has_deficit(imbalance)
deba@2440
   388
#endif
deba@2440
   389
    {
deba@2440
   390
      // Checking the sum of supply values
deba@2440
   391
      Supply sum = 0;
deba@2440
   392
      for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
deba@2440
   393
      valid_supply = sum == 0;
deba@2440
   394
    }
deba@2440
   395
deba@2440
   396
    /// \brief Simple constructor of the class (with lower bounds).
deba@2440
   397
    ///
deba@2440
   398
    /// Simple constructor of the class (with lower bounds).
deba@2440
   399
    ///
deba@2440
   400
    /// \param _graph The directed graph the algorithm runs on.
deba@2440
   401
    /// \param _lower The lower bounds of the edges.
deba@2440
   402
    /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440
   403
    /// \param _cost The cost (length) values of the edges.
deba@2440
   404
    /// \param _s The source node.
deba@2440
   405
    /// \param _t The target node.
deba@2440
   406
    /// \param _flow_value The required amount of flow from node \c _s
deba@2440
   407
    /// to node \c _t (i.e. the supply of \c _s and the demand of
deba@2440
   408
    /// \c _t).
deba@2440
   409
    CapacityScaling( const Graph &_graph,
deba@2440
   410
		     const LowerMap &_lower,
deba@2440
   411
		     const CapacityMap &_capacity,
deba@2440
   412
		     const CostMap &_cost,
deba@2440
   413
		     Node _s, Node _t,
deba@2440
   414
		     Supply _flow_value ) :
deba@2440
   415
      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
deba@2440
   416
      supply(_graph), flow(_graph, 0), potential(_graph, 0),
deba@2440
   417
      res_graph(_graph, capacity, flow),
deba@2440
   418
      red_cost(res_graph, cost, potential), imbalance(_graph),
deba@2440
   419
#ifdef WITH_SCALING
deba@2440
   420
      delta(0), delta_filter(res_graph, delta),
deba@2440
   421
      dres_graph(res_graph, delta_filter),
deba@2440
   422
      dijkstra(dres_graph, red_cost), pred(dres_graph),
deba@2440
   423
      delta_deficit(imbalance, delta)
deba@2440
   424
#else //WITHOUT_CAPACITY_SCALING
deba@2440
   425
      dijkstra(res_graph, red_cost), pred(res_graph),
deba@2440
   426
      has_deficit(imbalance)
deba@2440
   427
#endif
deba@2440
   428
    {
deba@2440
   429
      // Removing nonzero lower bounds
deba@2440
   430
      capacity = subMap(_capacity, _lower);
deba@2440
   431
      for (NodeIt n(graph); n != INVALID; ++n) {
deba@2440
   432
	Supply s = 0;
deba@2440
   433
	if (n == _s) s =  _flow_value;
deba@2440
   434
	if (n == _t) s = -_flow_value;
deba@2440
   435
	for (InEdgeIt e(graph, n); e != INVALID; ++e)
deba@2440
   436
	  s += _lower[e];
deba@2440
   437
	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
deba@2440
   438
	  s -= _lower[e];
deba@2440
   439
	supply[n] = imbalance[n] = s;
deba@2440
   440
      }
deba@2440
   441
      valid_supply = true;
deba@2440
   442
    }
deba@2440
   443
deba@2440
   444
    /// \brief Simple constructor of the class (without lower bounds).
deba@2440
   445
    ///
deba@2440
   446
    /// Simple constructor of the class (without lower bounds).
deba@2440
   447
    ///
deba@2440
   448
    /// \param _graph The directed graph the algorithm runs on.
deba@2440
   449
    /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440
   450
    /// \param _cost The cost (length) values of the edges.
deba@2440
   451
    /// \param _s The source node.
deba@2440
   452
    /// \param _t The target node.
deba@2440
   453
    /// \param _flow_value The required amount of flow from node \c _s
deba@2440
   454
    /// to node \c _t (i.e. the supply of \c _s and the demand of
deba@2440
   455
    /// \c _t).
deba@2440
   456
    CapacityScaling( const Graph &_graph,
deba@2440
   457
		     const CapacityMap &_capacity,
deba@2440
   458
		     const CostMap &_cost,
deba@2440
   459
		     Node _s, Node _t,
deba@2440
   460
		     Supply _flow_value ) :
deba@2440
   461
      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
deba@2440
   462
      supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
deba@2440
   463
      res_graph(_graph, capacity, flow),
deba@2440
   464
      red_cost(res_graph, cost, potential), imbalance(_graph),
deba@2440
   465
#ifdef WITH_SCALING
deba@2440
   466
      delta(0), delta_filter(res_graph, delta),
deba@2440
   467
      dres_graph(res_graph, delta_filter),
deba@2440
   468
      dijkstra(dres_graph, red_cost), pred(dres_graph),
deba@2440
   469
      delta_deficit(imbalance, delta)
deba@2440
   470
#else //WITHOUT_CAPACITY_SCALING
deba@2440
   471
      dijkstra(res_graph, red_cost), pred(res_graph),
deba@2440
   472
      has_deficit(imbalance)
deba@2440
   473
#endif
deba@2440
   474
    {
deba@2440
   475
      supply[_s] =  _flow_value;
deba@2440
   476
      supply[_t] = -_flow_value;
deba@2440
   477
      valid_supply = true;
deba@2440
   478
    }
deba@2440
   479
deba@2440
   480
    /// \brief Returns a const reference to the flow map.
deba@2440
   481
    ///
deba@2440
   482
    /// Returns a const reference to the flow map.
deba@2440
   483
    ///
deba@2440
   484
    /// \pre \ref run() must be called before using this function.
deba@2440
   485
    const FlowMap& flowMap() const {
deba@2440
   486
      return flow;
deba@2440
   487
    }
deba@2440
   488
deba@2440
   489
    /// \brief Returns a const reference to the potential map (the dual
deba@2440
   490
    /// solution).
deba@2440
   491
    ///
deba@2440
   492
    /// Returns a const reference to the potential map (the dual
deba@2440
   493
    /// solution).
deba@2440
   494
    ///
deba@2440
   495
    /// \pre \ref run() must be called before using this function.
deba@2440
   496
    const PotentialMap& potentialMap() const {
deba@2440
   497
      return potential;
deba@2440
   498
    }
deba@2440
   499
deba@2440
   500
    /// \brief Returns the total cost of the found flow.
deba@2440
   501
    ///
deba@2440
   502
    /// Returns the total cost of the found flow. The complexity of the
deba@2440
   503
    /// function is \f$ O(e) \f$.
deba@2440
   504
    ///
deba@2440
   505
    /// \pre \ref run() must be called before using this function.
deba@2440
   506
    Cost totalCost() const {
deba@2440
   507
      Cost c = 0;
deba@2440
   508
      for (EdgeIt e(graph); e != INVALID; ++e)
deba@2440
   509
	c += flow[e] * cost[e];
deba@2440
   510
      return c;
deba@2440
   511
    }
deba@2440
   512
deba@2440
   513
    /// \brief Runs the successive shortest path algorithm.
deba@2440
   514
    ///
deba@2440
   515
    /// Runs the successive shortest path algorithm.
deba@2440
   516
    ///
deba@2440
   517
    /// \return \c true if a feasible flow can be found.
deba@2440
   518
    bool run() {
deba@2440
   519
      return init() && start();
deba@2440
   520
    }
deba@2440
   521
deba@2440
   522
  protected:
deba@2440
   523
deba@2440
   524
    /// \brief Initializes the algorithm.
deba@2440
   525
    bool init() {
deba@2440
   526
      if (!valid_supply) return false;
deba@2440
   527
deba@2440
   528
      // Initalizing Dijkstra class
deba@2440
   529
      updater.potentialMap(potential);
deba@2440
   530
      dijkstra.distMap(updater).predMap(pred);
deba@2440
   531
deba@2440
   532
#ifdef WITH_SCALING
deba@2440
   533
      // Initilaizing delta value
deba@2440
   534
      Capacity max_cap = 0;
deba@2440
   535
      for (EdgeIt e(graph); e != INVALID; ++e) {
deba@2440
   536
	if (capacity[e] > max_cap) max_cap = capacity[e];
deba@2440
   537
      }
deba@2440
   538
      for (delta = 1; 2 * delta < max_cap; delta *= 2) ;
deba@2440
   539
#endif
deba@2440
   540
      return true;
deba@2440
   541
    }
deba@2440
   542
deba@2440
   543
#ifdef WITH_SCALING
deba@2440
   544
    /// \brief Executes the capacity scaling version of the successive
deba@2440
   545
    /// shortest path algorithm.
deba@2440
   546
    bool start() {
deba@2440
   547
      typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
deba@2440
   548
      typedef typename DeltaResGraph::Edge DeltaResEdge;
deba@2440
   549
deba@2440
   550
      // Processing capacity scaling phases
deba@2440
   551
      ResNode s, t;
deba@2440
   552
      for ( ; delta >= 1; delta = delta < 4 && delta > 1 ?
deba@2440
   553
				  1 : delta / 4 )
deba@2440
   554
      {
deba@2440
   555
	// Saturating edges not satisfying the optimality condition
deba@2440
   556
	Capacity r;
deba@2440
   557
	for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {
deba@2440
   558
	  if (red_cost[e] < 0) {
deba@2440
   559
	    res_graph.augment(e, r = res_graph.rescap(e));
deba@2440
   560
	    imbalance[dres_graph.target(e)] += r;
deba@2440
   561
	    imbalance[dres_graph.source(e)] -= r;
deba@2440
   562
	  }
deba@2440
   563
	}
deba@2440
   564
deba@2440
   565
	// Finding excess nodes
deba@2440
   566
	excess_nodes.clear();
deba@2440
   567
	for (ResNodeIt n(res_graph); n != INVALID; ++n) {
deba@2440
   568
	  if (imbalance[n] >= delta) excess_nodes.push_back(n);
deba@2440
   569
	}
deba@2440
   570
	next_node = 0;
deba@2440
   571
deba@2440
   572
	// Finding successive shortest paths
deba@2440
   573
	while (next_node < excess_nodes.size()) {
deba@2440
   574
	  // Running Dijkstra
deba@2440
   575
	  s = excess_nodes[next_node];
deba@2440
   576
	  updater.init();
deba@2440
   577
	  dijkstra.init();
deba@2440
   578
	  dijkstra.addSource(s);
deba@2440
   579
	  if ((t = dijkstra.start(delta_deficit)) == INVALID) {
deba@2440
   580
	    if (delta > 1) {
deba@2440
   581
	      ++next_node;
deba@2440
   582
	      continue;
deba@2440
   583
	    }
deba@2440
   584
	    return false;
deba@2440
   585
	  }
deba@2440
   586
deba@2440
   587
	  // Updating node potentials
deba@2440
   588
	  updater.update();
deba@2440
   589
deba@2440
   590
	  // Augment along a shortest path from s to t
deba@2440
   591
	  Capacity d = imbalance[s] < -imbalance[t] ?
deba@2440
   592
	    imbalance[s] : -imbalance[t];
deba@2440
   593
	  ResNode u = t;
deba@2440
   594
	  ResEdge e;
deba@2440
   595
	  if (d > delta) {
deba@2440
   596
	    while ((e = pred[u]) != INVALID) {
deba@2440
   597
	      if (res_graph.rescap(e) < d) d = res_graph.rescap(e);
deba@2440
   598
	      u = dres_graph.source(e);
deba@2440
   599
	    }
deba@2440
   600
	  }
deba@2440
   601
	  u = t;
deba@2440
   602
	  while ((e = pred[u]) != INVALID) {
deba@2440
   603
	    res_graph.augment(e, d);
deba@2440
   604
	    u = dres_graph.source(e);
deba@2440
   605
	  }
deba@2440
   606
	  imbalance[s] -= d;
deba@2440
   607
	  imbalance[t] += d;
deba@2440
   608
	  if (imbalance[s] < delta) ++next_node;
deba@2440
   609
	}
deba@2440
   610
      }
deba@2440
   611
deba@2440
   612
      // Handling nonzero lower bounds
deba@2440
   613
      if (lower) {
deba@2440
   614
	for (EdgeIt e(graph); e != INVALID; ++e)
deba@2440
   615
	  flow[e] += (*lower)[e];
deba@2440
   616
      }
deba@2440
   617
      return true;
deba@2440
   618
    }
deba@2440
   619
deba@2440
   620
#else //WITHOUT_CAPACITY_SCALING
deba@2440
   621
    /// \brief Executes the successive shortest path algorithm without
deba@2440
   622
    /// capacity scaling.
deba@2440
   623
    bool start() {
deba@2440
   624
      // Finding excess nodes
deba@2440
   625
      for (ResNodeIt n(res_graph); n != INVALID; ++n) {
deba@2440
   626
	if (imbalance[n] > 0) excess_nodes.push_back(n);
deba@2440
   627
      }
deba@2440
   628
      if (excess_nodes.size() == 0) return true;
deba@2440
   629
      next_node = 0;
deba@2440
   630
deba@2440
   631
      // Finding successive shortest paths
deba@2440
   632
      ResNode s, t;
deba@2440
   633
      while ( imbalance[excess_nodes[next_node]] > 0 ||
deba@2440
   634
	      ++next_node < excess_nodes.size() )
deba@2440
   635
      {
deba@2440
   636
	// Running Dijkstra
deba@2440
   637
	s = excess_nodes[next_node];
deba@2440
   638
	updater.init();
deba@2440
   639
	dijkstra.init();
deba@2440
   640
	dijkstra.addSource(s);
deba@2440
   641
	if ((t = dijkstra.start(has_deficit)) == INVALID)
deba@2440
   642
	  return false;
deba@2440
   643
deba@2440
   644
	// Updating node potentials
deba@2440
   645
	updater.update();
deba@2440
   646
deba@2440
   647
	// Augmenting along a shortest path from s to t
deba@2440
   648
	Capacity delta = imbalance[s] < -imbalance[t] ?
deba@2440
   649
	  imbalance[s] : -imbalance[t];
deba@2440
   650
	ResNode u = t;
deba@2440
   651
	ResEdge e;
deba@2440
   652
	while ((e = pred[u]) != INVALID) {
deba@2440
   653
	  if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);
deba@2440
   654
	  u = res_graph.source(e);
deba@2440
   655
	}
deba@2440
   656
	u = t;
deba@2440
   657
	while ((e = pred[u]) != INVALID) {
deba@2440
   658
	  res_graph.augment(e, delta);
deba@2440
   659
	  u = res_graph.source(e);
deba@2440
   660
	}
deba@2440
   661
	imbalance[s] -= delta;
deba@2440
   662
	imbalance[t] += delta;
deba@2440
   663
      }
deba@2440
   664
deba@2440
   665
      // Handling nonzero lower bounds
deba@2440
   666
      if (lower) {
deba@2440
   667
	for (EdgeIt e(graph); e != INVALID; ++e)
deba@2440
   668
	  flow[e] += (*lower)[e];
deba@2440
   669
      }
deba@2440
   670
      return true;
deba@2440
   671
    }
deba@2440
   672
#endif
deba@2440
   673
deba@2440
   674
  }; //class CapacityScaling
deba@2440
   675
deba@2440
   676
  ///@}
deba@2440
   677
deba@2440
   678
} //namespace lemon
deba@2440
   679
deba@2440
   680
#endif //LEMON_CAPACITY_SCALING_H