Small bug fixes.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2007
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_CAPACITY_SCALING_H
20 #define LEMON_CAPACITY_SCALING_H
22 /// \ingroup min_cost_flow
25 /// \brief The capacity scaling algorithm for finding a minimum cost
29 #include <lemon/dijkstra.h>
30 #include <lemon/graph_adaptor.h>
34 //#define _DEBUG_ITER_
38 /// \addtogroup min_cost_flow
41 /// \brief Implementation of the capacity scaling version of the
42 /// successive shortest path algorithm for finding a minimum cost flow.
44 /// \ref lemon::CapacityScaling "CapacityScaling" implements a
45 /// capacity scaling algorithm for finding a minimum cost flow.
47 /// \param Graph The directed graph type the algorithm runs on.
48 /// \param LowerMap The type of the lower bound map.
49 /// \param CapacityMap The type of the capacity (upper bound) map.
50 /// \param CostMap The type of the cost (length) map.
51 /// \param SupplyMap The type of the supply map.
54 /// - Edge capacities and costs should be nonnegative integers.
55 /// However \c CostMap::Value should be signed type.
56 /// - Supply values should be integers.
57 /// - \c LowerMap::Value must be convertible to
58 /// \c CapacityMap::Value and \c CapacityMap::Value must be
59 /// convertible to \c SupplyMap::Value.
61 /// \author Peter Kovacs
63 template < typename Graph,
64 typename LowerMap = typename Graph::template EdgeMap<int>,
65 typename CapacityMap = LowerMap,
66 typename CostMap = typename Graph::template EdgeMap<int>,
67 typename SupplyMap = typename Graph::template NodeMap
68 <typename CapacityMap::Value> >
71 typedef typename Graph::Node Node;
72 typedef typename Graph::NodeIt NodeIt;
73 typedef typename Graph::Edge Edge;
74 typedef typename Graph::EdgeIt EdgeIt;
75 typedef typename Graph::InEdgeIt InEdgeIt;
76 typedef typename Graph::OutEdgeIt OutEdgeIt;
78 typedef typename LowerMap::Value Lower;
79 typedef typename CapacityMap::Value Capacity;
80 typedef typename CostMap::Value Cost;
81 typedef typename SupplyMap::Value Supply;
82 typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
83 typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
85 typedef ResGraphAdaptor< const Graph, Capacity,
86 CapacityRefMap, CapacityRefMap > ResGraph;
87 typedef typename ResGraph::Node ResNode;
88 typedef typename ResGraph::NodeIt ResNodeIt;
89 typedef typename ResGraph::Edge ResEdge;
90 typedef typename ResGraph::EdgeIt ResEdgeIt;
94 /// \brief The type of the flow map.
95 typedef CapacityRefMap FlowMap;
96 /// \brief The type of the potential map.
97 typedef typename Graph::template NodeMap<Cost> PotentialMap;
101 /// \brief Map adaptor class for handling reduced edge costs.
102 class ReducedCostMap : public MapBase<ResEdge, Cost>
107 const CostMap &cost_map;
108 const PotentialMap &pot_map;
112 typedef typename MapBase<ResEdge, Cost>::Value Value;
113 typedef typename MapBase<ResEdge, Cost>::Key Key;
115 ReducedCostMap( const ResGraph &_gr,
116 const CostMap &_cost,
117 const PotentialMap &_pot ) :
118 gr(_gr), cost_map(_cost), pot_map(_pot) {}
120 Value operator[](const Key &e) const {
121 return ResGraph::forward(e) ?
122 cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :
123 -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
126 }; //class ReducedCostMap
128 /// \brief Map adaptor for \ref lemon::Dijkstra "Dijkstra" class to
129 /// update node potentials.
130 class PotentialUpdateMap : public MapBase<ResNode, Cost>
135 typedef std::pair<ResNode, Cost> Pair;
136 std::vector<Pair> data;
140 typedef typename MapBase<ResNode, Cost>::Value Value;
141 typedef typename MapBase<ResNode, Cost>::Key Key;
143 void potentialMap(PotentialMap &_pot) {
151 void set(const Key &n, const Value &v) {
152 data.push_back(Pair(n, v));
156 Cost val = data[data.size()-1].second;
157 for (int i = 0; i < data.size()-1; ++i)
158 (*pot)[data[i].first] += val - data[i].second;
161 }; //class PotentialUpdateMap
164 /// \brief Map adaptor class for identifing deficit nodes.
165 class DeficitBoolMap : public MapBase<ResNode, bool>
169 const SupplyRefMap &imb;
170 const Capacity δ
174 DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :
175 imb(_imb), delta(_delta) {}
177 bool operator[](const ResNode &n) const {
178 return imb[n] <= -delta;
181 }; //class DeficitBoolMap
183 /// \brief Map adaptor class for filtering edges with at least
184 /// \c delta residual capacity
185 class DeltaFilterMap : public MapBase<ResEdge, bool>
190 const Capacity δ
194 typedef typename MapBase<ResEdge, Cost>::Value Value;
195 typedef typename MapBase<ResEdge, Cost>::Key Key;
197 DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :
198 gr(_gr), delta(_delta) {}
200 Value operator[](const Key &e) const {
201 return gr.rescap(e) >= delta;
204 }; //class DeltaFilterMap
206 typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>
209 /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
210 class ResDijkstraTraits :
211 public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>
215 typedef PotentialUpdateMap DistMap;
217 static DistMap *createDistMap(const DeltaResGraph&) {
218 return new DistMap();
221 }; //class ResDijkstraTraits
223 #else //WITHOUT_CAPACITY_SCALING
224 /// \brief Map adaptor class for identifing deficit nodes.
225 class DeficitBoolMap : public MapBase<ResNode, bool>
229 const SupplyRefMap &imb;
233 DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}
235 bool operator[](const ResNode &n) const {
239 }; //class DeficitBoolMap
241 /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
242 class ResDijkstraTraits :
243 public DijkstraDefaultTraits<ResGraph, ReducedCostMap>
247 typedef PotentialUpdateMap DistMap;
249 static DistMap *createDistMap(const ResGraph&) {
250 return new DistMap();
253 }; //class ResDijkstraTraits
258 /// \brief The directed graph the algorithm runs on.
260 /// \brief The original lower bound map.
261 const LowerMap *lower;
262 /// \brief The modified capacity map.
263 CapacityRefMap capacity;
264 /// \brief The cost map.
266 /// \brief The modified supply map.
268 /// \brief The sum of supply values equals zero.
271 /// \brief The edge map of the current flow.
273 /// \brief The potential node map.
274 PotentialMap potential;
275 /// \brief The residual graph.
277 /// \brief The reduced cost map.
278 ReducedCostMap red_cost;
280 /// \brief The imbalance map.
281 SupplyRefMap imbalance;
282 /// \brief The excess nodes.
283 std::vector<ResNode> excess_nodes;
284 /// \brief The index of the next excess node.
288 typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>
290 /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
291 /// shortest paths in the residual graph with respect to the
292 /// reduced edge costs.
293 ResDijkstra dijkstra;
295 /// \brief The delta parameter used for capacity scaling.
297 /// \brief Edge filter map.
298 DeltaFilterMap delta_filter;
299 /// \brief The delta residual graph.
300 DeltaResGraph dres_graph;
301 /// \brief Map for identifing deficit nodes.
302 DeficitBoolMap delta_deficit;
304 /// \brief The deficit nodes.
305 std::vector<ResNode> deficit_nodes;
307 #else //WITHOUT_CAPACITY_SCALING
308 typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
310 /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
311 /// shortest paths in the residual graph with respect to the
312 /// reduced edge costs.
313 ResDijkstra dijkstra;
314 /// \brief Map for identifing deficit nodes.
315 DeficitBoolMap has_deficit;
318 /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.
319 typename ResDijkstra::PredMap pred;
320 /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to
321 /// update node potentials.
322 PotentialUpdateMap updater;
326 /// \brief General constructor of the class (with lower bounds).
328 /// General constructor of the class (with lower bounds).
330 /// \param _graph The directed graph the algorithm runs on.
331 /// \param _lower The lower bounds of the edges.
332 /// \param _capacity The capacities (upper bounds) of the edges.
333 /// \param _cost The cost (length) values of the edges.
334 /// \param _supply The supply values of the nodes (signed).
335 CapacityScaling( const Graph &_graph,
336 const LowerMap &_lower,
337 const CapacityMap &_capacity,
338 const CostMap &_cost,
339 const SupplyMap &_supply ) :
340 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
341 supply(_graph), flow(_graph, 0), potential(_graph, 0),
342 res_graph(_graph, capacity, flow),
343 red_cost(res_graph, cost, potential), imbalance(_graph),
345 delta(0), delta_filter(res_graph, delta),
346 dres_graph(res_graph, delta_filter),
347 dijkstra(dres_graph, red_cost), pred(dres_graph),
348 delta_deficit(imbalance, delta)
349 #else //WITHOUT_CAPACITY_SCALING
350 dijkstra(res_graph, red_cost), pred(res_graph),
351 has_deficit(imbalance)
354 // Removing nonzero lower bounds
355 capacity = subMap(_capacity, _lower);
357 for (NodeIt n(graph); n != INVALID; ++n) {
358 Supply s = _supply[n];
359 for (InEdgeIt e(graph, n); e != INVALID; ++e)
361 for (OutEdgeIt e(graph, n); e != INVALID; ++e)
363 supply[n] = imbalance[n] = s;
366 valid_supply = sum == 0;
369 /// \brief General constructor of the class (without lower bounds).
371 /// General constructor of the class (without lower bounds).
373 /// \param _graph The directed graph the algorithm runs on.
374 /// \param _capacity The capacities (upper bounds) of the edges.
375 /// \param _cost The cost (length) values of the edges.
376 /// \param _supply The supply values of the nodes (signed).
377 CapacityScaling( const Graph &_graph,
378 const CapacityMap &_capacity,
379 const CostMap &_cost,
380 const SupplyMap &_supply ) :
381 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
382 supply(_supply), flow(_graph, 0), potential(_graph, 0),
383 res_graph(_graph, capacity, flow),
384 red_cost(res_graph, cost, potential), imbalance(_graph),
386 delta(0), delta_filter(res_graph, delta),
387 dres_graph(res_graph, delta_filter),
388 dijkstra(dres_graph, red_cost), pred(dres_graph),
389 delta_deficit(imbalance, delta)
390 #else //WITHOUT_CAPACITY_SCALING
391 dijkstra(res_graph, red_cost), pred(res_graph),
392 has_deficit(imbalance)
395 // Checking the sum of supply values
397 for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
398 valid_supply = sum == 0;
401 /// \brief Simple constructor of the class (with lower bounds).
403 /// Simple constructor of the class (with lower bounds).
405 /// \param _graph The directed graph the algorithm runs on.
406 /// \param _lower The lower bounds of the edges.
407 /// \param _capacity The capacities (upper bounds) of the edges.
408 /// \param _cost The cost (length) values of the edges.
409 /// \param _s The source node.
410 /// \param _t The target node.
411 /// \param _flow_value The required amount of flow from node \c _s
412 /// to node \c _t (i.e. the supply of \c _s and the demand of
414 CapacityScaling( const Graph &_graph,
415 const LowerMap &_lower,
416 const CapacityMap &_capacity,
417 const CostMap &_cost,
419 Supply _flow_value ) :
420 graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
421 supply(_graph), flow(_graph, 0), potential(_graph, 0),
422 res_graph(_graph, capacity, flow),
423 red_cost(res_graph, cost, potential), imbalance(_graph),
425 delta(0), delta_filter(res_graph, delta),
426 dres_graph(res_graph, delta_filter),
427 dijkstra(dres_graph, red_cost), pred(dres_graph),
428 delta_deficit(imbalance, delta)
429 #else //WITHOUT_CAPACITY_SCALING
430 dijkstra(res_graph, red_cost), pred(res_graph),
431 has_deficit(imbalance)
434 // Removing nonzero lower bounds
435 capacity = subMap(_capacity, _lower);
436 for (NodeIt n(graph); n != INVALID; ++n) {
438 if (n == _s) s = _flow_value;
439 if (n == _t) s = -_flow_value;
440 for (InEdgeIt e(graph, n); e != INVALID; ++e)
442 for (OutEdgeIt e(graph, n); e != INVALID; ++e)
444 supply[n] = imbalance[n] = s;
449 /// \brief Simple constructor of the class (without lower bounds).
451 /// Simple constructor of the class (without lower bounds).
453 /// \param _graph The directed graph the algorithm runs on.
454 /// \param _capacity The capacities (upper bounds) of the edges.
455 /// \param _cost The cost (length) values of the edges.
456 /// \param _s The source node.
457 /// \param _t The target node.
458 /// \param _flow_value The required amount of flow from node \c _s
459 /// to node \c _t (i.e. the supply of \c _s and the demand of
461 CapacityScaling( const Graph &_graph,
462 const CapacityMap &_capacity,
463 const CostMap &_cost,
465 Supply _flow_value ) :
466 graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
467 supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
468 res_graph(_graph, capacity, flow),
469 red_cost(res_graph, cost, potential), imbalance(_graph),
471 delta(0), delta_filter(res_graph, delta),
472 dres_graph(res_graph, delta_filter),
473 dijkstra(dres_graph, red_cost), pred(dres_graph),
474 delta_deficit(imbalance, delta)
475 #else //WITHOUT_CAPACITY_SCALING
476 dijkstra(res_graph, red_cost), pred(res_graph),
477 has_deficit(imbalance)
480 supply[_s] = _flow_value;
481 supply[_t] = -_flow_value;
485 /// \brief Returns a const reference to the flow map.
487 /// Returns a const reference to the flow map.
489 /// \pre \ref run() must be called before using this function.
490 const FlowMap& flowMap() const {
494 /// \brief Returns a const reference to the potential map (the dual
497 /// Returns a const reference to the potential map (the dual
500 /// \pre \ref run() must be called before using this function.
501 const PotentialMap& potentialMap() const {
505 /// \brief Returns the total cost of the found flow.
507 /// Returns the total cost of the found flow. The complexity of the
508 /// function is \f$ O(e) \f$.
510 /// \pre \ref run() must be called before using this function.
511 Cost totalCost() const {
513 for (EdgeIt e(graph); e != INVALID; ++e)
514 c += flow[e] * cost[e];
518 /// \brief Runs the algorithm.
520 /// Runs the algorithm.
522 /// \return \c true if a feasible flow can be found.
524 return init() && start();
529 /// \brief Initializes the algorithm.
531 if (!valid_supply) return false;
533 // Initalizing Dijkstra class
534 updater.potentialMap(potential);
535 dijkstra.distMap(updater).predMap(pred);
538 // Initilaizing delta value
539 Supply max_sup = 0, max_dem = 0;
540 for (NodeIt n(graph); n != INVALID; ++n) {
541 if (supply[n] > max_sup) max_sup = supply[n];
542 if (supply[n] < -max_dem) max_dem = -supply[n];
544 if (max_dem < max_sup) max_sup = max_dem;
545 for (delta = 1; 2 * delta < max_sup; delta *= 2) ;
551 /// \brief Executes the capacity scaling version of the successive
552 /// shortest path algorithm.
554 typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
555 typedef typename DeltaResGraph::Edge DeltaResEdge;
560 // Processing capacity scaling phases
562 for ( ; delta >= 1; delta = delta < 4 && delta > 1 ?
565 // Saturating edges not satisfying the optimality condition
567 for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {
568 if (red_cost[e] < 0) {
569 res_graph.augment(e, r = res_graph.rescap(e));
570 imbalance[dres_graph.target(e)] += r;
571 imbalance[dres_graph.source(e)] -= r;
575 // Finding excess nodes
576 excess_nodes.clear();
577 deficit_nodes.clear();
578 for (ResNodeIt n(res_graph); n != INVALID; ++n) {
579 if (imbalance[n] >= delta) excess_nodes.push_back(n);
580 if (imbalance[n] <= -delta) deficit_nodes.push_back(n);
584 // Finding shortest paths
585 while (next_node < excess_nodes.size()) {
586 // Checking deficit nodes
588 bool delta_def = false;
589 for (int i = 0; i < deficit_nodes.size(); ++i) {
590 if (imbalance[deficit_nodes[i]] <= -delta) {
595 if (!delta_def) break;
599 s = excess_nodes[next_node];
602 dijkstra.addSource(s);
606 if ((t = dijkstra.start(delta_deficit)) == INVALID) {
614 // Updating node potentials
617 // Augment along a shortest path from s to t
618 Capacity d = imbalance[s] < -imbalance[t] ?
619 imbalance[s] : -imbalance[t];
623 while ((e = pred[u]) != INVALID) {
624 if (res_graph.rescap(e) < d) d = res_graph.rescap(e);
625 u = dres_graph.source(e);
629 while ((e = pred[u]) != INVALID) {
630 res_graph.augment(e, d);
631 u = dres_graph.source(e);
635 if (imbalance[s] < delta) ++next_node;
639 std::cout << "Cost Scaling algorithm finished with running Dijkstra algorithm "
640 << dijk_num << " times." << std::endl;
643 // Handling nonzero lower bounds
645 for (EdgeIt e(graph); e != INVALID; ++e)
646 flow[e] += (*lower)[e];
651 #else //WITHOUT_CAPACITY_SCALING
652 /// \brief Executes the successive shortest path algorithm without
653 /// capacity scaling.
655 // Finding excess nodes
656 for (ResNodeIt n(res_graph); n != INVALID; ++n) {
657 if (imbalance[n] > 0) excess_nodes.push_back(n);
659 if (excess_nodes.size() == 0) return true;
662 // Finding shortest paths
664 while ( imbalance[excess_nodes[next_node]] > 0 ||
665 ++next_node < excess_nodes.size() )
668 s = excess_nodes[next_node];
671 dijkstra.addSource(s);
672 if ((t = dijkstra.start(has_deficit)) == INVALID)
675 // Updating node potentials
678 // Augmenting along a shortest path from s to t
679 Capacity delta = imbalance[s] < -imbalance[t] ?
680 imbalance[s] : -imbalance[t];
683 while ((e = pred[u]) != INVALID) {
684 if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);
685 u = res_graph.source(e);
688 while ((e = pred[u]) != INVALID) {
689 res_graph.augment(e, delta);
690 u = res_graph.source(e);
692 imbalance[s] -= delta;
693 imbalance[t] += delta;
696 // Handling nonzero lower bounds
698 for (EdgeIt e(graph); e != INVALID; ++e)
699 flow[e] += (*lower)[e];
705 }; //class CapacityScaling
711 #endif //LEMON_CAPACITY_SCALING_H