Redesigned CapacityScaling algorithm with almost the same interface.
The new version does not use the ResidualGraphAdaptor for performance reasons.
Scaling can be enabled and disabled with a parameter of the run() function.
1.1 --- a/lemon/capacity_scaling.h Wed Dec 05 12:57:24 2007 +0000
1.2 +++ b/lemon/capacity_scaling.h Wed Dec 05 13:03:19 2007 +0000
1.3 @@ -25,17 +25,9 @@
1.4 /// \brief The capacity scaling algorithm for finding a minimum cost
1.5 /// flow.
1.6
1.7 +#include <lemon/graph_adaptor.h>
1.8 +#include <lemon/bin_heap.h>
1.9 #include <vector>
1.10 -#include <lemon/graph_adaptor.h>
1.11 -#include <lemon/dijkstra.h>
1.12 -
1.13 -#define WITH_SCALING
1.14 -
1.15 -#ifdef WITH_SCALING
1.16 -#define SCALING_FACTOR 2
1.17 -#endif
1.18 -
1.19 -//#define _DEBUG_ITER_
1.20
1.21 namespace lemon {
1.22
1.23 @@ -46,9 +38,9 @@
1.24 /// successive shortest path algorithm for finding a minimum cost
1.25 /// flow.
1.26 ///
1.27 - /// \ref lemon::CapacityScaling "CapacityScaling" implements the
1.28 - /// capacity scaling version of the successive shortest path
1.29 - /// algorithm for finding a minimum cost flow.
1.30 + /// \ref CapacityScaling implements the capacity scaling version
1.31 + /// of the successive shortest path algorithm for finding a minimum
1.32 + /// cost flow.
1.33 ///
1.34 /// \param Graph The directed graph type the algorithm runs on.
1.35 /// \param LowerMap The type of the lower bound map.
1.36 @@ -58,20 +50,20 @@
1.37 ///
1.38 /// \warning
1.39 /// - Edge capacities and costs should be nonnegative integers.
1.40 - /// However \c CostMap::Value should be signed type.
1.41 + /// However \c CostMap::Value should be signed type.
1.42 /// - Supply values should be signed integers.
1.43 /// - \c LowerMap::Value must be convertible to
1.44 - /// \c CapacityMap::Value and \c CapacityMap::Value must be
1.45 - /// convertible to \c SupplyMap::Value.
1.46 + /// \c CapacityMap::Value and \c CapacityMap::Value must be
1.47 + /// convertible to \c SupplyMap::Value.
1.48 ///
1.49 /// \author Peter Kovacs
1.50
1.51 template < typename Graph,
1.52 - typename LowerMap = typename Graph::template EdgeMap<int>,
1.53 - typename CapacityMap = LowerMap,
1.54 - typename CostMap = typename Graph::template EdgeMap<int>,
1.55 - typename SupplyMap = typename Graph::template NodeMap
1.56 - <typename CapacityMap::Value> >
1.57 + typename LowerMap = typename Graph::template EdgeMap<int>,
1.58 + typename CapacityMap = LowerMap,
1.59 + typename CostMap = typename Graph::template EdgeMap<int>,
1.60 + typename SupplyMap = typename Graph::template NodeMap
1.61 + <typename CapacityMap::Value> >
1.62 class CapacityScaling
1.63 {
1.64 typedef typename Graph::Node Node;
1.65 @@ -87,16 +79,16 @@
1.66 typedef typename SupplyMap::Value Supply;
1.67 typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
1.68 typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
1.69 -
1.70 - typedef ResGraphAdaptor< const Graph, Capacity,
1.71 - CapacityRefMap, CapacityRefMap > ResGraph;
1.72 - typedef typename ResGraph::Node ResNode;
1.73 - typedef typename ResGraph::NodeIt ResNodeIt;
1.74 - typedef typename ResGraph::Edge ResEdge;
1.75 - typedef typename ResGraph::EdgeIt ResEdgeIt;
1.76 + typedef typename Graph::template NodeMap<Edge> PredMap;
1.77
1.78 public:
1.79
1.80 + /// \brief Type to enable or disable capacity scaling.
1.81 + enum ScalingEnum {
1.82 + WITH_SCALING = 0,
1.83 + WITHOUT_SCALING = -1
1.84 + };
1.85 +
1.86 /// \brief The type of the flow map.
1.87 typedef CapacityRefMap FlowMap;
1.88 /// \brief The type of the potential map.
1.89 @@ -104,151 +96,129 @@
1.90
1.91 protected:
1.92
1.93 - /// \brief Map adaptor class for handling reduced edge costs.
1.94 - class ReducedCostMap : public MapBase<ResEdge, Cost>
1.95 + /// \brief Special implementation of the \ref Dijkstra algorithm
1.96 + /// for finding shortest paths in the residual network of the graph
1.97 + /// with respect to the reduced edge costs and modifying the
1.98 + /// node potentials according to the distance of the nodes.
1.99 + class ResidualDijkstra
1.100 {
1.101 - private:
1.102 + typedef typename Graph::template NodeMap<Cost> CostNodeMap;
1.103 + typedef typename Graph::template NodeMap<Edge> PredMap;
1.104
1.105 - const ResGraph &gr;
1.106 - const CostMap &cost_map;
1.107 - const PotentialMap &pot_map;
1.108 + typedef typename Graph::template NodeMap<int> HeapCrossRef;
1.109 + typedef BinHeap<Cost, HeapCrossRef> Heap;
1.110 +
1.111 + protected:
1.112 +
1.113 + /// \brief The directed graph the algorithm runs on.
1.114 + const Graph &graph;
1.115 +
1.116 + /// \brief The flow map.
1.117 + const FlowMap &flow;
1.118 + /// \brief The residual capacity map.
1.119 + const CapacityRefMap &res_cap;
1.120 + /// \brief The cost map.
1.121 + const CostMap &cost;
1.122 + /// \brief The excess map.
1.123 + const SupplyRefMap &excess;
1.124 +
1.125 + /// \brief The potential map.
1.126 + PotentialMap &potential;
1.127 +
1.128 + /// \brief The distance map.
1.129 + CostNodeMap dist;
1.130 + /// \brief The map of predecessors edges.
1.131 + PredMap &pred;
1.132 + /// \brief The processed (i.e. permanently labeled) nodes.
1.133 + std::vector<Node> proc_nodes;
1.134
1.135 public:
1.136
1.137 - ReducedCostMap( const ResGraph &_gr,
1.138 - const CostMap &_cost,
1.139 - const PotentialMap &_pot ) :
1.140 - gr(_gr), cost_map(_cost), pot_map(_pot) {}
1.141 + /// \brief The constructor of the class.
1.142 + ResidualDijkstra( const Graph &_graph,
1.143 + const FlowMap &_flow,
1.144 + const CapacityRefMap &_res_cap,
1.145 + const CostMap &_cost,
1.146 + const SupplyMap &_excess,
1.147 + PotentialMap &_potential,
1.148 + PredMap &_pred ) :
1.149 + graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost),
1.150 + excess(_excess), potential(_potential), dist(_graph),
1.151 + pred(_pred)
1.152 + {}
1.153
1.154 - Cost operator[](const ResEdge &e) const {
1.155 - return ResGraph::forward(e) ?
1.156 - cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :
1.157 - -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
1.158 + /// \brief Runs the algorithm from the given source node.
1.159 + Node run(Node s, Capacity delta) {
1.160 + HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP);
1.161 + Heap heap(heap_cross_ref);
1.162 + heap.push(s, 0);
1.163 + pred[s] = INVALID;
1.164 + proc_nodes.clear();
1.165 +
1.166 + // Processing nodes
1.167 + while (!heap.empty() && excess[heap.top()] > -delta) {
1.168 + Node u = heap.top(), v;
1.169 + Cost d = heap.prio() - potential[u], nd;
1.170 + dist[u] = heap.prio();
1.171 + heap.pop();
1.172 + proc_nodes.push_back(u);
1.173 +
1.174 + // Traversing outgoing edges
1.175 + for (OutEdgeIt e(graph, u); e != INVALID; ++e) {
1.176 + if (res_cap[e] >= delta) {
1.177 + v = graph.target(e);
1.178 + switch(heap.state(v)) {
1.179 + case Heap::PRE_HEAP:
1.180 + heap.push(v, d + cost[e] + potential[v]);
1.181 + pred[v] = e;
1.182 + break;
1.183 + case Heap::IN_HEAP:
1.184 + nd = d + cost[e] + potential[v];
1.185 + if (nd < heap[v]) {
1.186 + heap.decrease(v, nd);
1.187 + pred[v] = e;
1.188 + }
1.189 + break;
1.190 + case Heap::POST_HEAP:
1.191 + break;
1.192 + }
1.193 + }
1.194 + }
1.195 +
1.196 + // Traversing incoming edges
1.197 + for (InEdgeIt e(graph, u); e != INVALID; ++e) {
1.198 + if (flow[e] >= delta) {
1.199 + v = graph.source(e);
1.200 + switch(heap.state(v)) {
1.201 + case Heap::PRE_HEAP:
1.202 + heap.push(v, d - cost[e] + potential[v]);
1.203 + pred[v] = e;
1.204 + break;
1.205 + case Heap::IN_HEAP:
1.206 + nd = d - cost[e] + potential[v];
1.207 + if (nd < heap[v]) {
1.208 + heap.decrease(v, nd);
1.209 + pred[v] = e;
1.210 + }
1.211 + break;
1.212 + case Heap::POST_HEAP:
1.213 + break;
1.214 + }
1.215 + }
1.216 + }
1.217 + }
1.218 + if (heap.empty()) return INVALID;
1.219 +
1.220 + // Updating potentials of processed nodes
1.221 + Node t = heap.top();
1.222 + Cost dt = heap.prio();
1.223 + for (int i = 0; i < proc_nodes.size(); ++i)
1.224 + potential[proc_nodes[i]] -= dist[proc_nodes[i]] - dt;
1.225 +
1.226 + return t;
1.227 }
1.228
1.229 - }; //class ReducedCostMap
1.230 -
1.231 - /// \brief Map class for the \ref lemon::Dijkstra "Dijkstra"
1.232 - /// algorithm to update node potentials.
1.233 - class PotentialUpdateMap : public MapBase<ResNode, Cost>
1.234 - {
1.235 - private:
1.236 -
1.237 - PotentialMap *pot;
1.238 - typedef std::pair<ResNode, Cost> Pair;
1.239 - std::vector<Pair> data;
1.240 -
1.241 - public:
1.242 -
1.243 - void potentialMap(PotentialMap &_pot) {
1.244 - pot = &_pot;
1.245 - }
1.246 -
1.247 - void init() {
1.248 - data.clear();
1.249 - }
1.250 -
1.251 - void set(const ResNode &n, const Cost &v) {
1.252 - data.push_back(Pair(n, v));
1.253 - }
1.254 -
1.255 - void update() {
1.256 - Cost val = data[data.size()-1].second;
1.257 - for (int i = 0; i < data.size()-1; ++i)
1.258 - (*pot)[data[i].first] += val - data[i].second;
1.259 - }
1.260 -
1.261 - }; //class PotentialUpdateMap
1.262 -
1.263 -#ifdef WITH_SCALING
1.264 - /// \brief Map adaptor class for identifing deficit nodes.
1.265 - class DeficitBoolMap : public MapBase<ResNode, bool>
1.266 - {
1.267 - private:
1.268 -
1.269 - const SupplyRefMap &imb;
1.270 - const Capacity δ
1.271 -
1.272 - public:
1.273 -
1.274 - DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :
1.275 - imb(_imb), delta(_delta) {}
1.276 -
1.277 - bool operator[](const ResNode &n) const {
1.278 - return imb[n] <= -delta;
1.279 - }
1.280 -
1.281 - }; //class DeficitBoolMap
1.282 -
1.283 - /// \brief Map adaptor class for filtering edges with at least
1.284 - /// \c delta residual capacity.
1.285 - class DeltaFilterMap : public MapBase<ResEdge, bool>
1.286 - {
1.287 - private:
1.288 -
1.289 - const ResGraph &gr;
1.290 - const Capacity δ
1.291 -
1.292 - public:
1.293 -
1.294 - DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :
1.295 - gr(_gr), delta(_delta) {}
1.296 -
1.297 - bool operator[](const ResEdge &e) const {
1.298 - return gr.rescap(e) >= delta;
1.299 - }
1.300 -
1.301 - }; //class DeltaFilterMap
1.302 -
1.303 - typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>
1.304 - DeltaResGraph;
1.305 -
1.306 - /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
1.307 - class ResDijkstraTraits :
1.308 - public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>
1.309 - {
1.310 - public:
1.311 -
1.312 - typedef PotentialUpdateMap DistMap;
1.313 -
1.314 - static DistMap *createDistMap(const DeltaResGraph&) {
1.315 - return new DistMap();
1.316 - }
1.317 -
1.318 - }; //class ResDijkstraTraits
1.319 -
1.320 -#else //WITHOUT_CAPACITY_SCALING
1.321 - /// \brief Map adaptor class for identifing deficit nodes.
1.322 - class DeficitBoolMap : public MapBase<ResNode, bool>
1.323 - {
1.324 - private:
1.325 -
1.326 - const SupplyRefMap &imb;
1.327 -
1.328 - public:
1.329 -
1.330 - DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}
1.331 -
1.332 - bool operator[](const ResNode &n) const {
1.333 - return imb[n] < 0;
1.334 - }
1.335 -
1.336 - }; //class DeficitBoolMap
1.337 -
1.338 - /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
1.339 - class ResDijkstraTraits :
1.340 - public DijkstraDefaultTraits<ResGraph, ReducedCostMap>
1.341 - {
1.342 - public:
1.343 -
1.344 - typedef PotentialUpdateMap DistMap;
1.345 -
1.346 - static DistMap *createDistMap(const ResGraph&) {
1.347 - return new DistMap();
1.348 - }
1.349 -
1.350 - }; //class ResDijkstraTraits
1.351 -#endif
1.352 + }; //class ResidualDijkstra
1.353
1.354 protected:
1.355
1.356 @@ -269,57 +239,30 @@
1.357 FlowMap flow;
1.358 /// \brief The potential node map.
1.359 PotentialMap potential;
1.360 - /// \brief The residual graph.
1.361 - ResGraph res_graph;
1.362 - /// \brief The reduced cost map.
1.363 - ReducedCostMap red_cost;
1.364
1.365 - /// \brief The imbalance map.
1.366 - SupplyRefMap imbalance;
1.367 - /// \brief The excess nodes.
1.368 - std::vector<ResNode> excess_nodes;
1.369 + /// \brief The residual capacity map.
1.370 + CapacityRefMap res_cap;
1.371 + /// \brief The excess map.
1.372 + SupplyRefMap excess;
1.373 + /// \brief The excess nodes (i.e. nodes with positive excess).
1.374 + std::vector<Node> excess_nodes;
1.375 /// \brief The index of the next excess node.
1.376 int next_node;
1.377
1.378 -#ifdef WITH_SCALING
1.379 - typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>
1.380 - ResDijkstra;
1.381 - /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
1.382 - /// shortest paths in the residual graph with respect to the
1.383 - /// reduced edge costs.
1.384 - ResDijkstra dijkstra;
1.385 -
1.386 + /// \brief The scaling status (enabled or disabled).
1.387 + ScalingEnum scaling;
1.388 /// \brief The delta parameter used for capacity scaling.
1.389 Capacity delta;
1.390 - /// \brief Edge filter map for filtering edges with at least
1.391 - /// \c delta residual capacity.
1.392 - DeltaFilterMap delta_filter;
1.393 - /// \brief The delta residual graph (i.e. the subgraph of the
1.394 - /// residual graph consisting of edges with at least \c delta
1.395 - /// residual capacity).
1.396 - DeltaResGraph dres_graph;
1.397 - /// \brief Map for identifing deficit nodes.
1.398 - DeficitBoolMap delta_deficit;
1.399 + /// \brief The maximum number of phases.
1.400 + Capacity phase_num;
1.401 + /// \brief The deficit nodes.
1.402 + std::vector<Node> deficit_nodes;
1.403
1.404 - /// \brief The deficit nodes.
1.405 - std::vector<ResNode> deficit_nodes;
1.406 -
1.407 -#else //WITHOUT_CAPACITY_SCALING
1.408 - typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
1.409 - ResDijkstra;
1.410 - /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
1.411 - /// shortest paths in the residual graph with respect to the
1.412 - /// reduced edge costs.
1.413 - ResDijkstra dijkstra;
1.414 - /// \brief Map for identifing deficit nodes.
1.415 - DeficitBoolMap has_deficit;
1.416 -#endif
1.417 -
1.418 - /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.
1.419 - typename ResDijkstra::PredMap pred;
1.420 - /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to
1.421 - /// update node potentials.
1.422 - PotentialUpdateMap updater;
1.423 + /// \brief Implementation of the \ref Dijkstra algorithm for
1.424 + /// finding augmenting shortest paths in the residual network.
1.425 + ResidualDijkstra dijkstra;
1.426 + /// \brief The map of predecessors edges.
1.427 + PredMap pred;
1.428
1.429 public :
1.430
1.431 @@ -333,35 +276,27 @@
1.432 /// \param _cost The cost (length) values of the edges.
1.433 /// \param _supply The supply values of the nodes (signed).
1.434 CapacityScaling( const Graph &_graph,
1.435 - const LowerMap &_lower,
1.436 - const CapacityMap &_capacity,
1.437 - const CostMap &_cost,
1.438 - const SupplyMap &_supply ) :
1.439 + const LowerMap &_lower,
1.440 + const CapacityMap &_capacity,
1.441 + const CostMap &_cost,
1.442 + const SupplyMap &_supply ) :
1.443 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
1.444 supply(_graph), flow(_graph, 0), potential(_graph, 0),
1.445 - res_graph(_graph, capacity, flow),
1.446 - red_cost(res_graph, cost, potential), imbalance(_graph),
1.447 -#ifdef WITH_SCALING
1.448 - delta(0), delta_filter(res_graph, delta),
1.449 - dres_graph(res_graph, delta_filter),
1.450 - dijkstra(dres_graph, red_cost), pred(dres_graph),
1.451 - delta_deficit(imbalance, delta)
1.452 -#else //WITHOUT_CAPACITY_SCALING
1.453 - dijkstra(res_graph, red_cost), pred(res_graph),
1.454 - has_deficit(imbalance)
1.455 -#endif
1.456 + res_cap(_graph), excess(_graph), pred(_graph),
1.457 + dijkstra(graph, flow, res_cap, cost, excess, potential, pred)
1.458 {
1.459 // Removing nonzero lower bounds
1.460 capacity = subMap(_capacity, _lower);
1.461 + res_cap = capacity;
1.462 Supply sum = 0;
1.463 for (NodeIt n(graph); n != INVALID; ++n) {
1.464 - Supply s = _supply[n];
1.465 - for (InEdgeIt e(graph, n); e != INVALID; ++e)
1.466 - s += _lower[e];
1.467 - for (OutEdgeIt e(graph, n); e != INVALID; ++e)
1.468 - s -= _lower[e];
1.469 - supply[n] = s;
1.470 - sum += s;
1.471 + Supply s = _supply[n];
1.472 + for (InEdgeIt e(graph, n); e != INVALID; ++e)
1.473 + s += _lower[e];
1.474 + for (OutEdgeIt e(graph, n); e != INVALID; ++e)
1.475 + s -= _lower[e];
1.476 + supply[n] = s;
1.477 + sum += s;
1.478 }
1.479 valid_supply = sum == 0;
1.480 }
1.481 @@ -375,22 +310,13 @@
1.482 /// \param _cost The cost (length) values of the edges.
1.483 /// \param _supply The supply values of the nodes (signed).
1.484 CapacityScaling( const Graph &_graph,
1.485 - const CapacityMap &_capacity,
1.486 - const CostMap &_cost,
1.487 - const SupplyMap &_supply ) :
1.488 + const CapacityMap &_capacity,
1.489 + const CostMap &_cost,
1.490 + const SupplyMap &_supply ) :
1.491 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
1.492 supply(_supply), flow(_graph, 0), potential(_graph, 0),
1.493 - res_graph(_graph, capacity, flow),
1.494 - red_cost(res_graph, cost, potential), imbalance(_graph),
1.495 -#ifdef WITH_SCALING
1.496 - delta(0), delta_filter(res_graph, delta),
1.497 - dres_graph(res_graph, delta_filter),
1.498 - dijkstra(dres_graph, red_cost), pred(dres_graph),
1.499 - delta_deficit(imbalance, delta)
1.500 -#else //WITHOUT_CAPACITY_SCALING
1.501 - dijkstra(res_graph, red_cost), pred(res_graph),
1.502 - has_deficit(imbalance)
1.503 -#endif
1.504 + res_cap(_capacity), excess(_graph), pred(_graph),
1.505 + dijkstra(graph, flow, res_cap, cost, excess, potential)
1.506 {
1.507 // Checking the sum of supply values
1.508 Supply sum = 0;
1.509 @@ -412,36 +338,28 @@
1.510 /// to node \c _t (i.e. the supply of \c _s and the demand of
1.511 /// \c _t).
1.512 CapacityScaling( const Graph &_graph,
1.513 - const LowerMap &_lower,
1.514 - const CapacityMap &_capacity,
1.515 - const CostMap &_cost,
1.516 - Node _s, Node _t,
1.517 - Supply _flow_value ) :
1.518 + const LowerMap &_lower,
1.519 + const CapacityMap &_capacity,
1.520 + const CostMap &_cost,
1.521 + Node _s, Node _t,
1.522 + Supply _flow_value ) :
1.523 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
1.524 supply(_graph), flow(_graph, 0), potential(_graph, 0),
1.525 - res_graph(_graph, capacity, flow),
1.526 - red_cost(res_graph, cost, potential), imbalance(_graph),
1.527 -#ifdef WITH_SCALING
1.528 - delta(0), delta_filter(res_graph, delta),
1.529 - dres_graph(res_graph, delta_filter),
1.530 - dijkstra(dres_graph, red_cost), pred(dres_graph),
1.531 - delta_deficit(imbalance, delta)
1.532 -#else //WITHOUT_CAPACITY_SCALING
1.533 - dijkstra(res_graph, red_cost), pred(res_graph),
1.534 - has_deficit(imbalance)
1.535 -#endif
1.536 + res_cap(_graph), excess(_graph), pred(_graph),
1.537 + dijkstra(graph, flow, res_cap, cost, excess, potential)
1.538 {
1.539 // Removing nonzero lower bounds
1.540 capacity = subMap(_capacity, _lower);
1.541 + res_cap = capacity;
1.542 for (NodeIt n(graph); n != INVALID; ++n) {
1.543 - Supply s = 0;
1.544 - if (n == _s) s = _flow_value;
1.545 - if (n == _t) s = -_flow_value;
1.546 - for (InEdgeIt e(graph, n); e != INVALID; ++e)
1.547 - s += _lower[e];
1.548 - for (OutEdgeIt e(graph, n); e != INVALID; ++e)
1.549 - s -= _lower[e];
1.550 - supply[n] = s;
1.551 + Supply s = 0;
1.552 + if (n == _s) s = _flow_value;
1.553 + if (n == _t) s = -_flow_value;
1.554 + for (InEdgeIt e(graph, n); e != INVALID; ++e)
1.555 + s += _lower[e];
1.556 + for (OutEdgeIt e(graph, n); e != INVALID; ++e)
1.557 + s -= _lower[e];
1.558 + supply[n] = s;
1.559 }
1.560 valid_supply = true;
1.561 }
1.562 @@ -459,23 +377,14 @@
1.563 /// to node \c _t (i.e. the supply of \c _s and the demand of
1.564 /// \c _t).
1.565 CapacityScaling( const Graph &_graph,
1.566 - const CapacityMap &_capacity,
1.567 - const CostMap &_cost,
1.568 - Node _s, Node _t,
1.569 - Supply _flow_value ) :
1.570 + const CapacityMap &_capacity,
1.571 + const CostMap &_cost,
1.572 + Node _s, Node _t,
1.573 + Supply _flow_value ) :
1.574 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
1.575 supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
1.576 - res_graph(_graph, capacity, flow),
1.577 - red_cost(res_graph, cost, potential), imbalance(_graph),
1.578 -#ifdef WITH_SCALING
1.579 - delta(0), delta_filter(res_graph, delta),
1.580 - dres_graph(res_graph, delta_filter),
1.581 - dijkstra(dres_graph, red_cost), pred(dres_graph),
1.582 - delta_deficit(imbalance, delta)
1.583 -#else //WITHOUT_CAPACITY_SCALING
1.584 - dijkstra(res_graph, red_cost), pred(res_graph),
1.585 - has_deficit(imbalance)
1.586 -#endif
1.587 + res_cap(_capacity), excess(_graph), pred(_graph),
1.588 + dijkstra(graph, flow, res_cap, cost, excess, potential)
1.589 {
1.590 supply[_s] = _flow_value;
1.591 supply[_t] = -_flow_value;
1.592 @@ -511,7 +420,7 @@
1.593 Cost totalCost() const {
1.594 Cost c = 0;
1.595 for (EdgeIt e(graph); e != INVALID; ++e)
1.596 - c += flow[e] * cost[e];
1.597 + c += flow[e] * cost[e];
1.598 return c;
1.599 }
1.600
1.601 @@ -519,190 +428,216 @@
1.602 ///
1.603 /// Runs the algorithm.
1.604 ///
1.605 + /// \param scaling_mode The scaling mode. In case of WITH_SCALING
1.606 + /// capacity scaling is enabled in the algorithm (this is the
1.607 + /// default value) otherwise it is disabled.
1.608 + /// If the maximum edge capacity and/or the amount of total supply
1.609 + /// is small, the algorithm could be faster without scaling.
1.610 + ///
1.611 /// \return \c true if a feasible flow can be found.
1.612 - bool run() {
1.613 - return init() && start();
1.614 + bool run(int scaling_mode = WITH_SCALING) {
1.615 + return init(scaling_mode) && start();
1.616 }
1.617
1.618 protected:
1.619
1.620 /// \brief Initializes the algorithm.
1.621 - bool init() {
1.622 + bool init(int scaling_mode) {
1.623 if (!valid_supply) return false;
1.624 - imbalance = supply;
1.625 + excess = supply;
1.626
1.627 - // Initalizing Dijkstra class
1.628 - updater.potentialMap(potential);
1.629 - dijkstra.distMap(updater).predMap(pred);
1.630 -
1.631 -#ifdef WITH_SCALING
1.632 // Initilaizing delta value
1.633 - Supply max_sup = 0, max_dem = 0;
1.634 - for (NodeIt n(graph); n != INVALID; ++n) {
1.635 - if (supply[n] > max_sup) max_sup = supply[n];
1.636 - if (supply[n] < -max_dem) max_dem = -supply[n];
1.637 - }
1.638 - if (max_dem < max_sup) max_sup = max_dem;
1.639 - for ( delta = 1; SCALING_FACTOR * delta < max_sup;
1.640 - delta *= SCALING_FACTOR ) ;
1.641 -#endif
1.642 - return true;
1.643 - }
1.644 -
1.645 -#ifdef WITH_SCALING
1.646 - /// \brief Executes the capacity scaling version of the successive
1.647 - /// shortest path algorithm.
1.648 - bool start() {
1.649 - typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
1.650 - typedef typename DeltaResGraph::Edge DeltaResEdge;
1.651 -#ifdef _DEBUG_ITER_
1.652 - int dijk_num = 0;
1.653 -#endif
1.654 -
1.655 - // Processing capacity scaling phases
1.656 - ResNode s, t;
1.657 - for ( ; delta >= 1; delta = delta < SCALING_FACTOR && delta > 1 ?
1.658 - 1 : delta / SCALING_FACTOR )
1.659 - {
1.660 - // Saturating edges not satisfying the optimality condition
1.661 - Capacity r;
1.662 - for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {
1.663 - if (red_cost[e] < 0) {
1.664 - res_graph.augment(e, r = res_graph.rescap(e));
1.665 - imbalance[dres_graph.source(e)] -= r;
1.666 - imbalance[dres_graph.target(e)] += r;
1.667 - }
1.668 - }
1.669 -
1.670 - // Finding excess nodes
1.671 - excess_nodes.clear();
1.672 - deficit_nodes.clear();
1.673 - for (ResNodeIt n(res_graph); n != INVALID; ++n) {
1.674 - if (imbalance[n] >= delta) excess_nodes.push_back(n);
1.675 - if (imbalance[n] <= -delta) deficit_nodes.push_back(n);
1.676 - }
1.677 - next_node = 0;
1.678 -
1.679 - // Finding shortest paths
1.680 - while (next_node < excess_nodes.size()) {
1.681 - // Checking deficit nodes
1.682 - if (delta > 1) {
1.683 - bool delta_def = false;
1.684 - for (int i = 0; i < deficit_nodes.size(); ++i) {
1.685 - if (imbalance[deficit_nodes[i]] <= -delta) {
1.686 - delta_def = true;
1.687 - break;
1.688 - }
1.689 - }
1.690 - if (!delta_def) break;
1.691 - }
1.692 -
1.693 - // Running Dijkstra
1.694 - s = excess_nodes[next_node];
1.695 - updater.init();
1.696 - dijkstra.init();
1.697 - dijkstra.addSource(s);
1.698 -#ifdef _DEBUG_ITER_
1.699 - ++dijk_num;
1.700 -#endif
1.701 - if ((t = dijkstra.start(delta_deficit)) == INVALID) {
1.702 - if (delta > 1) {
1.703 - ++next_node;
1.704 - continue;
1.705 - }
1.706 - return false;
1.707 - }
1.708 -
1.709 - // Updating node potentials
1.710 - updater.update();
1.711 -
1.712 - // Augment along a shortest path from s to t
1.713 - Capacity d = imbalance[s] < -imbalance[t] ?
1.714 - imbalance[s] : -imbalance[t];
1.715 - ResNode u = t;
1.716 - ResEdge e;
1.717 - if (d > delta) {
1.718 - while ((e = pred[u]) != INVALID) {
1.719 - if (res_graph.rescap(e) < d) d = res_graph.rescap(e);
1.720 - u = dres_graph.source(e);
1.721 - }
1.722 - }
1.723 - u = t;
1.724 - while ((e = pred[u]) != INVALID) {
1.725 - res_graph.augment(e, d);
1.726 - u = dres_graph.source(e);
1.727 - }
1.728 - imbalance[s] -= d;
1.729 - imbalance[t] += d;
1.730 - if (imbalance[s] < delta) ++next_node;
1.731 - }
1.732 - }
1.733 -#ifdef _DEBUG_ITER_
1.734 - std::cout << "Capacity Scaling algorithm finished with running Dijkstra algorithm "
1.735 - << dijk_num << " times." << std::endl;
1.736 -#endif
1.737 -
1.738 - // Handling nonzero lower bounds
1.739 - if (lower) {
1.740 - for (EdgeIt e(graph); e != INVALID; ++e)
1.741 - flow[e] += (*lower)[e];
1.742 + if (scaling_mode == WITH_SCALING) {
1.743 + // With scaling
1.744 + Supply max_sup = 0, max_dem = 0;
1.745 + for (NodeIt n(graph); n != INVALID; ++n) {
1.746 + if ( supply[n] > max_sup) max_sup = supply[n];
1.747 + if (-supply[n] > max_dem) max_dem = -supply[n];
1.748 + }
1.749 + if (max_dem < max_sup) max_sup = max_dem;
1.750 + phase_num = 0;
1.751 + for (delta = 1; 2 * delta <= max_sup; delta *= 2)
1.752 + ++phase_num;
1.753 + } else {
1.754 + // Without scaling
1.755 + delta = 1;
1.756 }
1.757 return true;
1.758 }
1.759
1.760 -#else //WITHOUT_CAPACITY_SCALING
1.761 + /// \brief Executes the algorithm.
1.762 + bool start() {
1.763 + if (delta > 1)
1.764 + return startWithScaling();
1.765 + else
1.766 + return startWithoutScaling();
1.767 + }
1.768 +
1.769 + /// \brief Executes the capacity scaling version of the successive
1.770 + /// shortest path algorithm.
1.771 + bool startWithScaling() {
1.772 + // Processing capacity scaling phases
1.773 + Node s, t;
1.774 + int phase_cnt = 0;
1.775 + int factor = 4;
1.776 + while (true) {
1.777 + // Saturating all edges not satisfying the optimality condition
1.778 + for (EdgeIt e(graph); e != INVALID; ++e) {
1.779 + Node u = graph.source(e), v = graph.target(e);
1.780 + Cost c = cost[e] - potential[u] + potential[v];
1.781 + if (c < 0 && res_cap[e] >= delta) {
1.782 + excess[u] -= res_cap[e];
1.783 + excess[v] += res_cap[e];
1.784 + flow[e] = capacity[e];
1.785 + res_cap[e] = 0;
1.786 + }
1.787 + else if (c > 0 && flow[e] >= delta) {
1.788 + excess[u] += flow[e];
1.789 + excess[v] -= flow[e];
1.790 + flow[e] = 0;
1.791 + res_cap[e] = capacity[e];
1.792 + }
1.793 + }
1.794 +
1.795 + // Finding excess nodes and deficit nodes
1.796 + excess_nodes.clear();
1.797 + deficit_nodes.clear();
1.798 + for (NodeIt n(graph); n != INVALID; ++n) {
1.799 + if (excess[n] >= delta) excess_nodes.push_back(n);
1.800 + if (excess[n] <= -delta) deficit_nodes.push_back(n);
1.801 + }
1.802 + next_node = 0;
1.803 +
1.804 + // Finding augmenting shortest paths
1.805 + while (next_node < excess_nodes.size()) {
1.806 + // Checking deficit nodes
1.807 + if (delta > 1) {
1.808 + bool delta_deficit = false;
1.809 + for (int i = 0; i < deficit_nodes.size(); ++i) {
1.810 + if (excess[deficit_nodes[i]] <= -delta) {
1.811 + delta_deficit = true;
1.812 + break;
1.813 + }
1.814 + }
1.815 + if (!delta_deficit) break;
1.816 + }
1.817 +
1.818 + // Running Dijkstra
1.819 + s = excess_nodes[next_node];
1.820 + if ((t = dijkstra.run(s, delta)) == INVALID) {
1.821 + if (delta > 1) {
1.822 + ++next_node;
1.823 + continue;
1.824 + }
1.825 + return false;
1.826 + }
1.827 +
1.828 + // Augmenting along a shortest path from s to t.
1.829 + Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
1.830 + Node u = t;
1.831 + Edge e;
1.832 + if (d > delta) {
1.833 + while ((e = pred[u]) != INVALID) {
1.834 + Capacity rc;
1.835 + if (u == graph.target(e)) {
1.836 + rc = res_cap[e];
1.837 + u = graph.source(e);
1.838 + } else {
1.839 + rc = flow[e];
1.840 + u = graph.target(e);
1.841 + }
1.842 + if (rc < d) d = rc;
1.843 + }
1.844 + }
1.845 + u = t;
1.846 + while ((e = pred[u]) != INVALID) {
1.847 + if (u == graph.target(e)) {
1.848 + flow[e] += d;
1.849 + res_cap[e] -= d;
1.850 + u = graph.source(e);
1.851 + } else {
1.852 + flow[e] -= d;
1.853 + res_cap[e] += d;
1.854 + u = graph.target(e);
1.855 + }
1.856 + }
1.857 + excess[s] -= d;
1.858 + excess[t] += d;
1.859 +
1.860 + if (excess[s] < delta) ++next_node;
1.861 + }
1.862 +
1.863 + if (delta == 1) break;
1.864 + if (++phase_cnt > phase_num / 4) factor = 2;
1.865 + delta = delta <= factor ? 1 : delta / factor;
1.866 + }
1.867 +
1.868 + // Handling nonzero lower bounds
1.869 + if (lower) {
1.870 + for (EdgeIt e(graph); e != INVALID; ++e)
1.871 + flow[e] += (*lower)[e];
1.872 + }
1.873 + return true;
1.874 + }
1.875 +
1.876 /// \brief Executes the successive shortest path algorithm without
1.877 /// capacity scaling.
1.878 - bool start() {
1.879 + bool startWithoutScaling() {
1.880 // Finding excess nodes
1.881 - for (ResNodeIt n(res_graph); n != INVALID; ++n) {
1.882 - if (imbalance[n] > 0) excess_nodes.push_back(n);
1.883 + for (NodeIt n(graph); n != INVALID; ++n) {
1.884 + if (excess[n] > 0) excess_nodes.push_back(n);
1.885 }
1.886 if (excess_nodes.size() == 0) return true;
1.887 next_node = 0;
1.888
1.889 // Finding shortest paths
1.890 - ResNode s, t;
1.891 - while ( imbalance[excess_nodes[next_node]] > 0 ||
1.892 - ++next_node < excess_nodes.size() )
1.893 + Node s, t;
1.894 + while ( excess[excess_nodes[next_node]] > 0 ||
1.895 + ++next_node < excess_nodes.size() )
1.896 {
1.897 - // Running Dijkstra
1.898 - s = excess_nodes[next_node];
1.899 - updater.init();
1.900 - dijkstra.init();
1.901 - dijkstra.addSource(s);
1.902 - if ((t = dijkstra.start(has_deficit)) == INVALID)
1.903 - return false;
1.904 + // Running Dijkstra
1.905 + s = excess_nodes[next_node];
1.906 + if ((t = dijkstra.run(s, 1)) == INVALID)
1.907 + return false;
1.908
1.909 - // Updating node potentials
1.910 - updater.update();
1.911 -
1.912 - // Augmenting along a shortest path from s to t
1.913 - Capacity delta = imbalance[s] < -imbalance[t] ?
1.914 - imbalance[s] : -imbalance[t];
1.915 - ResNode u = t;
1.916 - ResEdge e;
1.917 - while ((e = pred[u]) != INVALID) {
1.918 - if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);
1.919 - u = res_graph.source(e);
1.920 - }
1.921 - u = t;
1.922 - while ((e = pred[u]) != INVALID) {
1.923 - res_graph.augment(e, delta);
1.924 - u = res_graph.source(e);
1.925 - }
1.926 - imbalance[s] -= delta;
1.927 - imbalance[t] += delta;
1.928 + // Augmenting along a shortest path from s to t
1.929 + Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
1.930 + Node u = t;
1.931 + Edge e;
1.932 + while ((e = pred[u]) != INVALID) {
1.933 + Capacity rc;
1.934 + if (u == graph.target(e)) {
1.935 + rc = res_cap[e];
1.936 + u = graph.source(e);
1.937 + } else {
1.938 + rc = flow[e];
1.939 + u = graph.target(e);
1.940 + }
1.941 + if (rc < d) d = rc;
1.942 + }
1.943 + u = t;
1.944 + while ((e = pred[u]) != INVALID) {
1.945 + if (u == graph.target(e)) {
1.946 + flow[e] += d;
1.947 + res_cap[e] -= d;
1.948 + u = graph.source(e);
1.949 + } else {
1.950 + flow[e] -= d;
1.951 + res_cap[e] += d;
1.952 + u = graph.target(e);
1.953 + }
1.954 + }
1.955 + excess[s] -= d;
1.956 + excess[t] += d;
1.957 }
1.958
1.959 // Handling nonzero lower bounds
1.960 if (lower) {
1.961 - for (EdgeIt e(graph); e != INVALID; ++e)
1.962 - flow[e] += (*lower)[e];
1.963 + for (EdgeIt e(graph); e != INVALID; ++e)
1.964 + flow[e] += (*lower)[e];
1.965 }
1.966 return true;
1.967 }
1.968 -#endif
1.969
1.970 }; //class CapacityScaling
1.971