1.1 --- a/lemon/capacity_scaling.h Thu Nov 12 23:17:34 2009 +0100
1.2 +++ b/lemon/capacity_scaling.h Thu Nov 12 23:26:13 2009 +0100
1.3 @@ -19,396 +19,417 @@
1.4 #ifndef LEMON_CAPACITY_SCALING_H
1.5 #define LEMON_CAPACITY_SCALING_H
1.6
1.7 -/// \ingroup min_cost_flow
1.8 +/// \ingroup min_cost_flow_algs
1.9 ///
1.10 /// \file
1.11 -/// \brief Capacity scaling algorithm for finding a minimum cost flow.
1.12 +/// \brief Capacity Scaling algorithm for finding a minimum cost flow.
1.13
1.14 #include <vector>
1.15 +#include <limits>
1.16 +#include <lemon/core.h>
1.17 #include <lemon/bin_heap.h>
1.18
1.19 namespace lemon {
1.20
1.21 - /// \addtogroup min_cost_flow
1.22 + /// \addtogroup min_cost_flow_algs
1.23 /// @{
1.24
1.25 - /// \brief Implementation of the capacity scaling algorithm for
1.26 - /// finding a minimum cost flow.
1.27 + /// \brief Implementation of the Capacity Scaling algorithm for
1.28 + /// finding a \ref min_cost_flow "minimum cost flow".
1.29 ///
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 + /// of the successive shortest path algorithm for finding a
1.34 + /// \ref min_cost_flow "minimum cost flow". It is an efficient dual
1.35 + /// solution method.
1.36 ///
1.37 - /// \tparam Digraph The digraph type the algorithm runs on.
1.38 - /// \tparam LowerMap The type of the lower bound map.
1.39 - /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.40 - /// \tparam CostMap The type of the cost (length) map.
1.41 - /// \tparam SupplyMap The type of the supply map.
1.42 + /// Most of the parameters of the problem (except for the digraph)
1.43 + /// can be given using separate functions, and the algorithm can be
1.44 + /// executed using the \ref run() function. If some parameters are not
1.45 + /// specified, then default values will be used.
1.46 ///
1.47 - /// \warning
1.48 - /// - Arc capacities and costs should be \e non-negative \e integers.
1.49 - /// - Supply values should be \e signed \e integers.
1.50 - /// - The value types of the maps should be convertible to each other.
1.51 - /// - \c CostMap::Value must be signed type.
1.52 + /// \tparam GR The digraph type the algorithm runs on.
1.53 + /// \tparam V The value type used for flow amounts, capacity bounds
1.54 + /// and supply values in the algorithm. By default it is \c int.
1.55 + /// \tparam C The value type used for costs and potentials in the
1.56 + /// algorithm. By default it is the same as \c V.
1.57 ///
1.58 - /// \author Peter Kovacs
1.59 - template < typename Digraph,
1.60 - typename LowerMap = typename Digraph::template ArcMap<int>,
1.61 - typename CapacityMap = typename Digraph::template ArcMap<int>,
1.62 - typename CostMap = typename Digraph::template ArcMap<int>,
1.63 - typename SupplyMap = typename Digraph::template NodeMap<int> >
1.64 + /// \warning Both value types must be signed and all input data must
1.65 + /// be integer.
1.66 + /// \warning This algorithm does not support negative costs for such
1.67 + /// arcs that have infinite upper bound.
1.68 + template <typename GR, typename V = int, typename C = V>
1.69 class CapacityScaling
1.70 {
1.71 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.72 + public:
1.73
1.74 - typedef typename CapacityMap::Value Capacity;
1.75 - typedef typename CostMap::Value Cost;
1.76 - typedef typename SupplyMap::Value Supply;
1.77 - typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
1.78 - typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
1.79 - typedef typename Digraph::template NodeMap<Arc> PredMap;
1.80 + /// The type of the flow amounts, capacity bounds and supply values
1.81 + typedef V Value;
1.82 + /// The type of the arc costs
1.83 + typedef C Cost;
1.84
1.85 public:
1.86
1.87 - /// The type of the flow map.
1.88 - typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.89 - /// The type of the potential map.
1.90 - typedef typename Digraph::template NodeMap<Cost> PotentialMap;
1.91 + /// \brief Problem type constants for the \c run() function.
1.92 + ///
1.93 + /// Enum type containing the problem type constants that can be
1.94 + /// returned by the \ref run() function of the algorithm.
1.95 + enum ProblemType {
1.96 + /// The problem has no feasible solution (flow).
1.97 + INFEASIBLE,
1.98 + /// The problem has optimal solution (i.e. it is feasible and
1.99 + /// bounded), and the algorithm has found optimal flow and node
1.100 + /// potentials (primal and dual solutions).
1.101 + OPTIMAL,
1.102 + /// The digraph contains an arc of negative cost and infinite
1.103 + /// upper bound. It means that the objective function is unbounded
1.104 + /// on that arc, however note that it could actually be bounded
1.105 + /// over the feasible flows, but this algroithm cannot handle
1.106 + /// these cases.
1.107 + UNBOUNDED
1.108 + };
1.109 +
1.110 + private:
1.111 +
1.112 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.113 +
1.114 + typedef std::vector<Arc> ArcVector;
1.115 + typedef std::vector<Node> NodeVector;
1.116 + typedef std::vector<int> IntVector;
1.117 + typedef std::vector<bool> BoolVector;
1.118 + typedef std::vector<Value> ValueVector;
1.119 + typedef std::vector<Cost> CostVector;
1.120
1.121 private:
1.122
1.123 - /// \brief Special implementation of the \ref Dijkstra algorithm
1.124 - /// for finding shortest paths in the residual network.
1.125 + // Data related to the underlying digraph
1.126 + const GR &_graph;
1.127 + int _node_num;
1.128 + int _arc_num;
1.129 + int _res_arc_num;
1.130 + int _root;
1.131 +
1.132 + // Parameters of the problem
1.133 + bool _have_lower;
1.134 + Value _sum_supply;
1.135 +
1.136 + // Data structures for storing the digraph
1.137 + IntNodeMap _node_id;
1.138 + IntArcMap _arc_idf;
1.139 + IntArcMap _arc_idb;
1.140 + IntVector _first_out;
1.141 + BoolVector _forward;
1.142 + IntVector _source;
1.143 + IntVector _target;
1.144 + IntVector _reverse;
1.145 +
1.146 + // Node and arc data
1.147 + ValueVector _lower;
1.148 + ValueVector _upper;
1.149 + CostVector _cost;
1.150 + ValueVector _supply;
1.151 +
1.152 + ValueVector _res_cap;
1.153 + CostVector _pi;
1.154 + ValueVector _excess;
1.155 + IntVector _excess_nodes;
1.156 + IntVector _deficit_nodes;
1.157 +
1.158 + Value _delta;
1.159 + int _phase_num;
1.160 + IntVector _pred;
1.161 +
1.162 + public:
1.163 +
1.164 + /// \brief Constant for infinite upper bounds (capacities).
1.165 ///
1.166 - /// \ref ResidualDijkstra is a special implementation of the
1.167 - /// \ref Dijkstra algorithm for finding shortest paths in the
1.168 - /// residual network of the digraph with respect to the reduced arc
1.169 - /// costs and modifying the node potentials according to the
1.170 - /// distance of the nodes.
1.171 + /// Constant for infinite upper bounds (capacities).
1.172 + /// It is \c std::numeric_limits<Value>::infinity() if available,
1.173 + /// \c std::numeric_limits<Value>::max() otherwise.
1.174 + const Value INF;
1.175 +
1.176 + private:
1.177 +
1.178 + // Special implementation of the Dijkstra algorithm for finding
1.179 + // shortest paths in the residual network of the digraph with
1.180 + // respect to the reduced arc costs and modifying the node
1.181 + // potentials according to the found distance labels.
1.182 class ResidualDijkstra
1.183 {
1.184 - typedef typename Digraph::template NodeMap<int> HeapCrossRef;
1.185 + typedef RangeMap<int> HeapCrossRef;
1.186 typedef BinHeap<Cost, HeapCrossRef> Heap;
1.187
1.188 private:
1.189
1.190 - // The digraph the algorithm runs on
1.191 - const Digraph &_graph;
1.192 -
1.193 - // The main maps
1.194 - const FlowMap &_flow;
1.195 - const CapacityArcMap &_res_cap;
1.196 - const CostMap &_cost;
1.197 - const SupplyNodeMap &_excess;
1.198 - PotentialMap &_potential;
1.199 -
1.200 - // The distance map
1.201 - PotentialMap _dist;
1.202 - // The pred arc map
1.203 - PredMap &_pred;
1.204 - // The processed (i.e. permanently labeled) nodes
1.205 - std::vector<Node> _proc_nodes;
1.206 -
1.207 + int _node_num;
1.208 + const IntVector &_first_out;
1.209 + const IntVector &_target;
1.210 + const CostVector &_cost;
1.211 + const ValueVector &_res_cap;
1.212 + const ValueVector &_excess;
1.213 + CostVector &_pi;
1.214 + IntVector &_pred;
1.215 +
1.216 + IntVector _proc_nodes;
1.217 + CostVector _dist;
1.218 +
1.219 public:
1.220
1.221 - /// Constructor.
1.222 - ResidualDijkstra( const Digraph &digraph,
1.223 - const FlowMap &flow,
1.224 - const CapacityArcMap &res_cap,
1.225 - const CostMap &cost,
1.226 - const SupplyMap &excess,
1.227 - PotentialMap &potential,
1.228 - PredMap &pred ) :
1.229 - _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost),
1.230 - _excess(excess), _potential(potential), _dist(digraph),
1.231 - _pred(pred)
1.232 + ResidualDijkstra(CapacityScaling& cs) :
1.233 + _node_num(cs._node_num), _first_out(cs._first_out),
1.234 + _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap),
1.235 + _excess(cs._excess), _pi(cs._pi), _pred(cs._pred),
1.236 + _dist(cs._node_num)
1.237 {}
1.238
1.239 - /// Run the algorithm from the given source node.
1.240 - Node run(Node s, Capacity delta = 1) {
1.241 - HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
1.242 + int run(int s, Value delta = 1) {
1.243 + HeapCrossRef heap_cross_ref(_node_num, Heap::PRE_HEAP);
1.244 Heap heap(heap_cross_ref);
1.245 heap.push(s, 0);
1.246 - _pred[s] = INVALID;
1.247 + _pred[s] = -1;
1.248 _proc_nodes.clear();
1.249
1.250 - // Processing nodes
1.251 + // Process nodes
1.252 while (!heap.empty() && _excess[heap.top()] > -delta) {
1.253 - Node u = heap.top(), v;
1.254 - Cost d = heap.prio() + _potential[u], nd;
1.255 + int u = heap.top(), v;
1.256 + Cost d = heap.prio() + _pi[u], dn;
1.257 _dist[u] = heap.prio();
1.258 + _proc_nodes.push_back(u);
1.259 heap.pop();
1.260 - _proc_nodes.push_back(u);
1.261
1.262 - // Traversing outgoing arcs
1.263 - for (OutArcIt e(_graph, u); e != INVALID; ++e) {
1.264 - if (_res_cap[e] >= delta) {
1.265 - v = _graph.target(e);
1.266 - switch(heap.state(v)) {
1.267 + // Traverse outgoing residual arcs
1.268 + for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
1.269 + if (_res_cap[a] < delta) continue;
1.270 + v = _target[a];
1.271 + switch (heap.state(v)) {
1.272 case Heap::PRE_HEAP:
1.273 - heap.push(v, d + _cost[e] - _potential[v]);
1.274 - _pred[v] = e;
1.275 + heap.push(v, d + _cost[a] - _pi[v]);
1.276 + _pred[v] = a;
1.277 break;
1.278 case Heap::IN_HEAP:
1.279 - nd = d + _cost[e] - _potential[v];
1.280 - if (nd < heap[v]) {
1.281 - heap.decrease(v, nd);
1.282 - _pred[v] = e;
1.283 + dn = d + _cost[a] - _pi[v];
1.284 + if (dn < heap[v]) {
1.285 + heap.decrease(v, dn);
1.286 + _pred[v] = a;
1.287 }
1.288 break;
1.289 case Heap::POST_HEAP:
1.290 break;
1.291 - }
1.292 - }
1.293 - }
1.294 -
1.295 - // Traversing incoming arcs
1.296 - for (InArcIt e(_graph, u); e != INVALID; ++e) {
1.297 - if (_flow[e] >= delta) {
1.298 - v = _graph.source(e);
1.299 - switch(heap.state(v)) {
1.300 - case Heap::PRE_HEAP:
1.301 - heap.push(v, d - _cost[e] - _potential[v]);
1.302 - _pred[v] = e;
1.303 - break;
1.304 - case Heap::IN_HEAP:
1.305 - nd = d - _cost[e] - _potential[v];
1.306 - if (nd < heap[v]) {
1.307 - heap.decrease(v, nd);
1.308 - _pred[v] = e;
1.309 - }
1.310 - break;
1.311 - case Heap::POST_HEAP:
1.312 - break;
1.313 - }
1.314 }
1.315 }
1.316 }
1.317 - if (heap.empty()) return INVALID;
1.318 + if (heap.empty()) return -1;
1.319
1.320 - // Updating potentials of processed nodes
1.321 - Node t = heap.top();
1.322 - Cost t_dist = heap.prio();
1.323 - for (int i = 0; i < int(_proc_nodes.size()); ++i)
1.324 - _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
1.325 + // Update potentials of processed nodes
1.326 + int t = heap.top();
1.327 + Cost dt = heap.prio();
1.328 + for (int i = 0; i < int(_proc_nodes.size()); ++i) {
1.329 + _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
1.330 + }
1.331
1.332 return t;
1.333 }
1.334
1.335 }; //class ResidualDijkstra
1.336
1.337 - private:
1.338 -
1.339 - // The digraph the algorithm runs on
1.340 - const Digraph &_graph;
1.341 - // The original lower bound map
1.342 - const LowerMap *_lower;
1.343 - // The modified capacity map
1.344 - CapacityArcMap _capacity;
1.345 - // The original cost map
1.346 - const CostMap &_cost;
1.347 - // The modified supply map
1.348 - SupplyNodeMap _supply;
1.349 - bool _valid_supply;
1.350 -
1.351 - // Arc map of the current flow
1.352 - FlowMap *_flow;
1.353 - bool _local_flow;
1.354 - // Node map of the current potentials
1.355 - PotentialMap *_potential;
1.356 - bool _local_potential;
1.357 -
1.358 - // The residual capacity map
1.359 - CapacityArcMap _res_cap;
1.360 - // The excess map
1.361 - SupplyNodeMap _excess;
1.362 - // The excess nodes (i.e. nodes with positive excess)
1.363 - std::vector<Node> _excess_nodes;
1.364 - // The deficit nodes (i.e. nodes with negative excess)
1.365 - std::vector<Node> _deficit_nodes;
1.366 -
1.367 - // The delta parameter used for capacity scaling
1.368 - Capacity _delta;
1.369 - // The maximum number of phases
1.370 - int _phase_num;
1.371 -
1.372 - // The pred arc map
1.373 - PredMap _pred;
1.374 - // Implementation of the Dijkstra algorithm for finding augmenting
1.375 - // shortest paths in the residual network
1.376 - ResidualDijkstra *_dijkstra;
1.377 -
1.378 public:
1.379
1.380 - /// \brief General constructor (with lower bounds).
1.381 + /// \brief Constructor.
1.382 ///
1.383 - /// General constructor (with lower bounds).
1.384 + /// The constructor of the class.
1.385 ///
1.386 - /// \param digraph The digraph the algorithm runs on.
1.387 - /// \param lower The lower bounds of the arcs.
1.388 - /// \param capacity The capacities (upper bounds) of the arcs.
1.389 - /// \param cost The cost (length) values of the arcs.
1.390 - /// \param supply The supply values of the nodes (signed).
1.391 - CapacityScaling( const Digraph &digraph,
1.392 - const LowerMap &lower,
1.393 - const CapacityMap &capacity,
1.394 - const CostMap &cost,
1.395 - const SupplyMap &supply ) :
1.396 - _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
1.397 - _supply(digraph), _flow(NULL), _local_flow(false),
1.398 - _potential(NULL), _local_potential(false),
1.399 - _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL)
1.400 + /// \param graph The digraph the algorithm runs on.
1.401 + CapacityScaling(const GR& graph) :
1.402 + _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
1.403 + INF(std::numeric_limits<Value>::has_infinity ?
1.404 + std::numeric_limits<Value>::infinity() :
1.405 + std::numeric_limits<Value>::max())
1.406 {
1.407 - Supply sum = 0;
1.408 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.409 - _supply[n] = supply[n];
1.410 - _excess[n] = supply[n];
1.411 - sum += supply[n];
1.412 + // Check the value types
1.413 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.414 + "The flow type of CapacityScaling must be signed");
1.415 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.416 + "The cost type of CapacityScaling must be signed");
1.417 +
1.418 + // Resize vectors
1.419 + _node_num = countNodes(_graph);
1.420 + _arc_num = countArcs(_graph);
1.421 + _res_arc_num = 2 * (_arc_num + _node_num);
1.422 + _root = _node_num;
1.423 + ++_node_num;
1.424 +
1.425 + _first_out.resize(_node_num + 1);
1.426 + _forward.resize(_res_arc_num);
1.427 + _source.resize(_res_arc_num);
1.428 + _target.resize(_res_arc_num);
1.429 + _reverse.resize(_res_arc_num);
1.430 +
1.431 + _lower.resize(_res_arc_num);
1.432 + _upper.resize(_res_arc_num);
1.433 + _cost.resize(_res_arc_num);
1.434 + _supply.resize(_node_num);
1.435 +
1.436 + _res_cap.resize(_res_arc_num);
1.437 + _pi.resize(_node_num);
1.438 + _excess.resize(_node_num);
1.439 + _pred.resize(_node_num);
1.440 +
1.441 + // Copy the graph
1.442 + int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
1.443 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.444 + _node_id[n] = i;
1.445 }
1.446 - _valid_supply = sum == 0;
1.447 + i = 0;
1.448 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.449 + _first_out[i] = j;
1.450 + for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.451 + _arc_idf[a] = j;
1.452 + _forward[j] = true;
1.453 + _source[j] = i;
1.454 + _target[j] = _node_id[_graph.runningNode(a)];
1.455 + }
1.456 + for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.457 + _arc_idb[a] = j;
1.458 + _forward[j] = false;
1.459 + _source[j] = i;
1.460 + _target[j] = _node_id[_graph.runningNode(a)];
1.461 + }
1.462 + _forward[j] = false;
1.463 + _source[j] = i;
1.464 + _target[j] = _root;
1.465 + _reverse[j] = k;
1.466 + _forward[k] = true;
1.467 + _source[k] = _root;
1.468 + _target[k] = i;
1.469 + _reverse[k] = j;
1.470 + ++j; ++k;
1.471 + }
1.472 + _first_out[i] = j;
1.473 + _first_out[_node_num] = k;
1.474 for (ArcIt a(_graph); a != INVALID; ++a) {
1.475 - _capacity[a] = capacity[a];
1.476 - _res_cap[a] = capacity[a];
1.477 + int fi = _arc_idf[a];
1.478 + int bi = _arc_idb[a];
1.479 + _reverse[fi] = bi;
1.480 + _reverse[bi] = fi;
1.481 }
1.482 -
1.483 - // Remove non-zero lower bounds
1.484 - typename LowerMap::Value lcap;
1.485 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.486 - if ((lcap = lower[e]) != 0) {
1.487 - _capacity[e] -= lcap;
1.488 - _res_cap[e] -= lcap;
1.489 - _supply[_graph.source(e)] -= lcap;
1.490 - _supply[_graph.target(e)] += lcap;
1.491 - _excess[_graph.source(e)] -= lcap;
1.492 - _excess[_graph.target(e)] += lcap;
1.493 - }
1.494 - }
1.495 - }
1.496 -/*
1.497 - /// \brief General constructor (without lower bounds).
1.498 - ///
1.499 - /// General constructor (without lower bounds).
1.500 - ///
1.501 - /// \param digraph The digraph the algorithm runs on.
1.502 - /// \param capacity The capacities (upper bounds) of the arcs.
1.503 - /// \param cost The cost (length) values of the arcs.
1.504 - /// \param supply The supply values of the nodes (signed).
1.505 - CapacityScaling( const Digraph &digraph,
1.506 - const CapacityMap &capacity,
1.507 - const CostMap &cost,
1.508 - const SupplyMap &supply ) :
1.509 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.510 - _supply(supply), _flow(NULL), _local_flow(false),
1.511 - _potential(NULL), _local_potential(false),
1.512 - _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL)
1.513 - {
1.514 - // Check the sum of supply values
1.515 - Supply sum = 0;
1.516 - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.517 - _valid_supply = sum == 0;
1.518 +
1.519 + // Reset parameters
1.520 + reset();
1.521 }
1.522
1.523 - /// \brief Simple constructor (with lower bounds).
1.524 + /// \name Parameters
1.525 + /// The parameters of the algorithm can be specified using these
1.526 + /// functions.
1.527 +
1.528 + /// @{
1.529 +
1.530 + /// \brief Set the lower bounds on the arcs.
1.531 ///
1.532 - /// Simple constructor (with lower bounds).
1.533 + /// This function sets the lower bounds on the arcs.
1.534 + /// If it is not used before calling \ref run(), the lower bounds
1.535 + /// will be set to zero on all arcs.
1.536 ///
1.537 - /// \param digraph The digraph the algorithm runs on.
1.538 - /// \param lower The lower bounds of the arcs.
1.539 - /// \param capacity The capacities (upper bounds) of the arcs.
1.540 - /// \param cost The cost (length) values of the arcs.
1.541 - /// \param s The source node.
1.542 - /// \param t The target node.
1.543 - /// \param flow_value The required amount of flow from node \c s
1.544 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.545 - CapacityScaling( const Digraph &digraph,
1.546 - const LowerMap &lower,
1.547 - const CapacityMap &capacity,
1.548 - const CostMap &cost,
1.549 - Node s, Node t,
1.550 - Supply flow_value ) :
1.551 - _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
1.552 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.553 - _potential(NULL), _local_potential(false),
1.554 - _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
1.555 - {
1.556 - // Remove non-zero lower bounds
1.557 - _supply[s] = _excess[s] = flow_value;
1.558 - _supply[t] = _excess[t] = -flow_value;
1.559 - typename LowerMap::Value lcap;
1.560 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.561 - if ((lcap = lower[e]) != 0) {
1.562 - _capacity[e] -= lcap;
1.563 - _res_cap[e] -= lcap;
1.564 - _supply[_graph.source(e)] -= lcap;
1.565 - _supply[_graph.target(e)] += lcap;
1.566 - _excess[_graph.source(e)] -= lcap;
1.567 - _excess[_graph.target(e)] += lcap;
1.568 - }
1.569 + /// \param map An arc map storing the lower bounds.
1.570 + /// Its \c Value type must be convertible to the \c Value type
1.571 + /// of the algorithm.
1.572 + ///
1.573 + /// \return <tt>(*this)</tt>
1.574 + template <typename LowerMap>
1.575 + CapacityScaling& lowerMap(const LowerMap& map) {
1.576 + _have_lower = true;
1.577 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.578 + _lower[_arc_idf[a]] = map[a];
1.579 + _lower[_arc_idb[a]] = map[a];
1.580 }
1.581 - _valid_supply = true;
1.582 - }
1.583 -
1.584 - /// \brief Simple constructor (without lower bounds).
1.585 - ///
1.586 - /// Simple constructor (without lower bounds).
1.587 - ///
1.588 - /// \param digraph The digraph the algorithm runs on.
1.589 - /// \param capacity The capacities (upper bounds) of the arcs.
1.590 - /// \param cost The cost (length) values of the arcs.
1.591 - /// \param s The source node.
1.592 - /// \param t The target node.
1.593 - /// \param flow_value The required amount of flow from node \c s
1.594 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.595 - CapacityScaling( const Digraph &digraph,
1.596 - const CapacityMap &capacity,
1.597 - const CostMap &cost,
1.598 - Node s, Node t,
1.599 - Supply flow_value ) :
1.600 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.601 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.602 - _potential(NULL), _local_potential(false),
1.603 - _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
1.604 - {
1.605 - _supply[s] = _excess[s] = flow_value;
1.606 - _supply[t] = _excess[t] = -flow_value;
1.607 - _valid_supply = true;
1.608 - }
1.609 -*/
1.610 - /// Destructor.
1.611 - ~CapacityScaling() {
1.612 - if (_local_flow) delete _flow;
1.613 - if (_local_potential) delete _potential;
1.614 - delete _dijkstra;
1.615 - }
1.616 -
1.617 - /// \brief Set the flow map.
1.618 - ///
1.619 - /// Set the flow map.
1.620 - ///
1.621 - /// \return \c (*this)
1.622 - CapacityScaling& flowMap(FlowMap &map) {
1.623 - if (_local_flow) {
1.624 - delete _flow;
1.625 - _local_flow = false;
1.626 - }
1.627 - _flow = ↦
1.628 return *this;
1.629 }
1.630
1.631 - /// \brief Set the potential map.
1.632 + /// \brief Set the upper bounds (capacities) on the arcs.
1.633 ///
1.634 - /// Set the potential map.
1.635 + /// This function sets the upper bounds (capacities) on the arcs.
1.636 + /// If it is not used before calling \ref run(), the upper bounds
1.637 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.638 + /// unbounded from above on each arc).
1.639 ///
1.640 - /// \return \c (*this)
1.641 - CapacityScaling& potentialMap(PotentialMap &map) {
1.642 - if (_local_potential) {
1.643 - delete _potential;
1.644 - _local_potential = false;
1.645 + /// \param map An arc map storing the upper bounds.
1.646 + /// Its \c Value type must be convertible to the \c Value type
1.647 + /// of the algorithm.
1.648 + ///
1.649 + /// \return <tt>(*this)</tt>
1.650 + template<typename UpperMap>
1.651 + CapacityScaling& upperMap(const UpperMap& map) {
1.652 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.653 + _upper[_arc_idf[a]] = map[a];
1.654 }
1.655 - _potential = ↦
1.656 return *this;
1.657 }
1.658
1.659 + /// \brief Set the costs of the arcs.
1.660 + ///
1.661 + /// This function sets the costs of the arcs.
1.662 + /// If it is not used before calling \ref run(), the costs
1.663 + /// will be set to \c 1 on all arcs.
1.664 + ///
1.665 + /// \param map An arc map storing the costs.
1.666 + /// Its \c Value type must be convertible to the \c Cost type
1.667 + /// of the algorithm.
1.668 + ///
1.669 + /// \return <tt>(*this)</tt>
1.670 + template<typename CostMap>
1.671 + CapacityScaling& costMap(const CostMap& map) {
1.672 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.673 + _cost[_arc_idf[a]] = map[a];
1.674 + _cost[_arc_idb[a]] = -map[a];
1.675 + }
1.676 + return *this;
1.677 + }
1.678 +
1.679 + /// \brief Set the supply values of the nodes.
1.680 + ///
1.681 + /// This function sets the supply values of the nodes.
1.682 + /// If neither this function nor \ref stSupply() is used before
1.683 + /// calling \ref run(), the supply of each node will be set to zero.
1.684 + ///
1.685 + /// \param map A node map storing the supply values.
1.686 + /// Its \c Value type must be convertible to the \c Value type
1.687 + /// of the algorithm.
1.688 + ///
1.689 + /// \return <tt>(*this)</tt>
1.690 + template<typename SupplyMap>
1.691 + CapacityScaling& supplyMap(const SupplyMap& map) {
1.692 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.693 + _supply[_node_id[n]] = map[n];
1.694 + }
1.695 + return *this;
1.696 + }
1.697 +
1.698 + /// \brief Set single source and target nodes and a supply value.
1.699 + ///
1.700 + /// This function sets a single source node and a single target node
1.701 + /// and the required flow value.
1.702 + /// If neither this function nor \ref supplyMap() is used before
1.703 + /// calling \ref run(), the supply of each node will be set to zero.
1.704 + ///
1.705 + /// Using this function has the same effect as using \ref supplyMap()
1.706 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.707 + /// assigned to \c t and all other nodes have zero supply value.
1.708 + ///
1.709 + /// \param s The source node.
1.710 + /// \param t The target node.
1.711 + /// \param k The required amount of flow from node \c s to node \c t
1.712 + /// (i.e. the supply of \c s and the demand of \c t).
1.713 + ///
1.714 + /// \return <tt>(*this)</tt>
1.715 + CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
1.716 + for (int i = 0; i != _node_num; ++i) {
1.717 + _supply[i] = 0;
1.718 + }
1.719 + _supply[_node_id[s]] = k;
1.720 + _supply[_node_id[t]] = -k;
1.721 + return *this;
1.722 + }
1.723 +
1.724 + /// @}
1.725 +
1.726 /// \name Execution control
1.727
1.728 /// @{
1.729 @@ -416,15 +437,88 @@
1.730 /// \brief Run the algorithm.
1.731 ///
1.732 /// This function runs the algorithm.
1.733 + /// The paramters can be specified using functions \ref lowerMap(),
1.734 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
1.735 + /// For example,
1.736 + /// \code
1.737 + /// CapacityScaling<ListDigraph> cs(graph);
1.738 + /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.739 + /// .supplyMap(sup).run();
1.740 + /// \endcode
1.741 + ///
1.742 + /// This function can be called more than once. All the parameters
1.743 + /// that have been given are kept for the next call, unless
1.744 + /// \ref reset() is called, thus only the modified parameters
1.745 + /// have to be set again. See \ref reset() for examples.
1.746 + /// However the underlying digraph must not be modified after this
1.747 + /// class have been constructed, since it copies the digraph.
1.748 ///
1.749 /// \param scaling Enable or disable capacity scaling.
1.750 - /// If the maximum arc capacity and/or the amount of total supply
1.751 + /// If the maximum upper bound and/or the amount of total supply
1.752 /// is rather small, the algorithm could be slightly faster without
1.753 /// scaling.
1.754 ///
1.755 - /// \return \c true if a feasible flow can be found.
1.756 - bool run(bool scaling = true) {
1.757 - return init(scaling) && start();
1.758 + /// \return \c INFEASIBLE if no feasible flow exists,
1.759 + /// \n \c OPTIMAL if the problem has optimal solution
1.760 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.761 + /// optimal flow and node potentials (primal and dual solutions),
1.762 + /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
1.763 + /// and infinite upper bound. It means that the objective function
1.764 + /// is unbounded on that arc, however note that it could actually be
1.765 + /// bounded over the feasible flows, but this algroithm cannot handle
1.766 + /// these cases.
1.767 + ///
1.768 + /// \see ProblemType
1.769 + ProblemType run(bool scaling = true) {
1.770 + ProblemType pt = init(scaling);
1.771 + if (pt != OPTIMAL) return pt;
1.772 + return start();
1.773 + }
1.774 +
1.775 + /// \brief Reset all the parameters that have been given before.
1.776 + ///
1.777 + /// This function resets all the paramaters that have been given
1.778 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.779 + /// \ref costMap(), \ref supplyMap(), \ref stSupply().
1.780 + ///
1.781 + /// It is useful for multiple run() calls. If this function is not
1.782 + /// used, all the parameters given before are kept for the next
1.783 + /// \ref run() call.
1.784 + /// However the underlying digraph must not be modified after this
1.785 + /// class have been constructed, since it copies and extends the graph.
1.786 + ///
1.787 + /// For example,
1.788 + /// \code
1.789 + /// CapacityScaling<ListDigraph> cs(graph);
1.790 + ///
1.791 + /// // First run
1.792 + /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.793 + /// .supplyMap(sup).run();
1.794 + ///
1.795 + /// // Run again with modified cost map (reset() is not called,
1.796 + /// // so only the cost map have to be set again)
1.797 + /// cost[e] += 100;
1.798 + /// cs.costMap(cost).run();
1.799 + ///
1.800 + /// // Run again from scratch using reset()
1.801 + /// // (the lower bounds will be set to zero on all arcs)
1.802 + /// cs.reset();
1.803 + /// cs.upperMap(capacity).costMap(cost)
1.804 + /// .supplyMap(sup).run();
1.805 + /// \endcode
1.806 + ///
1.807 + /// \return <tt>(*this)</tt>
1.808 + CapacityScaling& reset() {
1.809 + for (int i = 0; i != _node_num; ++i) {
1.810 + _supply[i] = 0;
1.811 + }
1.812 + for (int j = 0; j != _res_arc_num; ++j) {
1.813 + _lower[j] = 0;
1.814 + _upper[j] = INF;
1.815 + _cost[j] = _forward[j] ? 1 : -1;
1.816 + }
1.817 + _have_lower = false;
1.818 + return *this;
1.819 }
1.820
1.821 /// @}
1.822 @@ -432,97 +526,186 @@
1.823 /// \name Query Functions
1.824 /// The results of the algorithm can be obtained using these
1.825 /// functions.\n
1.826 - /// \ref lemon::CapacityScaling::run() "run()" must be called before
1.827 - /// using them.
1.828 + /// The \ref run() function must be called before using them.
1.829
1.830 /// @{
1.831
1.832 - /// \brief Return a const reference to the arc map storing the
1.833 - /// found flow.
1.834 + /// \brief Return the total cost of the found flow.
1.835 ///
1.836 - /// Return a const reference to the arc map storing the found flow.
1.837 + /// This function returns the total cost of the found flow.
1.838 + /// Its complexity is O(e).
1.839 + ///
1.840 + /// \note The return type of the function can be specified as a
1.841 + /// template parameter. For example,
1.842 + /// \code
1.843 + /// cs.totalCost<double>();
1.844 + /// \endcode
1.845 + /// It is useful if the total cost cannot be stored in the \c Cost
1.846 + /// type of the algorithm, which is the default return type of the
1.847 + /// function.
1.848 ///
1.849 /// \pre \ref run() must be called before using this function.
1.850 - const FlowMap& flowMap() const {
1.851 - return *_flow;
1.852 + template <typename Number>
1.853 + Number totalCost() const {
1.854 + Number c = 0;
1.855 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.856 + int i = _arc_idb[a];
1.857 + c += static_cast<Number>(_res_cap[i]) *
1.858 + (-static_cast<Number>(_cost[i]));
1.859 + }
1.860 + return c;
1.861 }
1.862
1.863 - /// \brief Return a const reference to the node map storing the
1.864 - /// found potentials (the dual solution).
1.865 - ///
1.866 - /// Return a const reference to the node map storing the found
1.867 - /// potentials (the dual solution).
1.868 - ///
1.869 - /// \pre \ref run() must be called before using this function.
1.870 - const PotentialMap& potentialMap() const {
1.871 - return *_potential;
1.872 +#ifndef DOXYGEN
1.873 + Cost totalCost() const {
1.874 + return totalCost<Cost>();
1.875 }
1.876 +#endif
1.877
1.878 /// \brief Return the flow on the given arc.
1.879 ///
1.880 - /// Return the flow on the given arc.
1.881 + /// This function returns the flow on the given arc.
1.882 ///
1.883 /// \pre \ref run() must be called before using this function.
1.884 - Capacity flow(const Arc& arc) const {
1.885 - return (*_flow)[arc];
1.886 + Value flow(const Arc& a) const {
1.887 + return _res_cap[_arc_idb[a]];
1.888 }
1.889
1.890 - /// \brief Return the potential of the given node.
1.891 + /// \brief Return the flow map (the primal solution).
1.892 ///
1.893 - /// Return the potential of the given node.
1.894 + /// This function copies the flow value on each arc into the given
1.895 + /// map. The \c Value type of the algorithm must be convertible to
1.896 + /// the \c Value type of the map.
1.897 ///
1.898 /// \pre \ref run() must be called before using this function.
1.899 - Cost potential(const Node& node) const {
1.900 - return (*_potential)[node];
1.901 + template <typename FlowMap>
1.902 + void flowMap(FlowMap &map) const {
1.903 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.904 + map.set(a, _res_cap[_arc_idb[a]]);
1.905 + }
1.906 }
1.907
1.908 - /// \brief Return the total cost of the found flow.
1.909 + /// \brief Return the potential (dual value) of the given node.
1.910 ///
1.911 - /// Return the total cost of the found flow. The complexity of the
1.912 - /// function is \f$ O(e) \f$.
1.913 + /// This function returns the potential (dual value) of the
1.914 + /// given node.
1.915 ///
1.916 /// \pre \ref run() must be called before using this function.
1.917 - Cost totalCost() const {
1.918 - Cost c = 0;
1.919 - for (ArcIt e(_graph); e != INVALID; ++e)
1.920 - c += (*_flow)[e] * _cost[e];
1.921 - return c;
1.922 + Cost potential(const Node& n) const {
1.923 + return _pi[_node_id[n]];
1.924 + }
1.925 +
1.926 + /// \brief Return the potential map (the dual solution).
1.927 + ///
1.928 + /// This function copies the potential (dual value) of each node
1.929 + /// into the given map.
1.930 + /// The \c Cost type of the algorithm must be convertible to the
1.931 + /// \c Value type of the map.
1.932 + ///
1.933 + /// \pre \ref run() must be called before using this function.
1.934 + template <typename PotentialMap>
1.935 + void potentialMap(PotentialMap &map) const {
1.936 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.937 + map.set(n, _pi[_node_id[n]]);
1.938 + }
1.939 }
1.940
1.941 /// @}
1.942
1.943 private:
1.944
1.945 - /// Initialize the algorithm.
1.946 - bool init(bool scaling) {
1.947 - if (!_valid_supply) return false;
1.948 + // Initialize the algorithm
1.949 + ProblemType init(bool scaling) {
1.950 + if (_node_num == 0) return INFEASIBLE;
1.951
1.952 - // Initializing maps
1.953 - if (!_flow) {
1.954 - _flow = new FlowMap(_graph);
1.955 - _local_flow = true;
1.956 + // Check the sum of supply values
1.957 + _sum_supply = 0;
1.958 + for (int i = 0; i != _root; ++i) {
1.959 + _sum_supply += _supply[i];
1.960 }
1.961 - if (!_potential) {
1.962 - _potential = new PotentialMap(_graph);
1.963 - _local_potential = true;
1.964 + if (_sum_supply > 0) return INFEASIBLE;
1.965 +
1.966 + // Initialize maps
1.967 + for (int i = 0; i != _root; ++i) {
1.968 + _pi[i] = 0;
1.969 + _excess[i] = _supply[i];
1.970 }
1.971 - for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
1.972 - for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
1.973
1.974 - _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
1.975 - _excess, *_potential, _pred );
1.976 + // Remove non-zero lower bounds
1.977 + if (_have_lower) {
1.978 + for (int i = 0; i != _root; ++i) {
1.979 + for (int j = _first_out[i]; j != _first_out[i+1]; ++j) {
1.980 + if (_forward[j]) {
1.981 + Value c = _lower[j];
1.982 + if (c >= 0) {
1.983 + _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF;
1.984 + } else {
1.985 + _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF;
1.986 + }
1.987 + _excess[i] -= c;
1.988 + _excess[_target[j]] += c;
1.989 + } else {
1.990 + _res_cap[j] = 0;
1.991 + }
1.992 + }
1.993 + }
1.994 + } else {
1.995 + for (int j = 0; j != _res_arc_num; ++j) {
1.996 + _res_cap[j] = _forward[j] ? _upper[j] : 0;
1.997 + }
1.998 + }
1.999
1.1000 - // Initializing delta value
1.1001 + // Handle negative costs
1.1002 + for (int u = 0; u != _root; ++u) {
1.1003 + for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
1.1004 + Value rc = _res_cap[a];
1.1005 + if (_cost[a] < 0 && rc > 0) {
1.1006 + if (rc == INF) return UNBOUNDED;
1.1007 + _excess[u] -= rc;
1.1008 + _excess[_target[a]] += rc;
1.1009 + _res_cap[a] = 0;
1.1010 + _res_cap[_reverse[a]] += rc;
1.1011 + }
1.1012 + }
1.1013 + }
1.1014 +
1.1015 + // Handle GEQ supply type
1.1016 + if (_sum_supply < 0) {
1.1017 + _pi[_root] = 0;
1.1018 + _excess[_root] = -_sum_supply;
1.1019 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.1020 + int u = _target[a];
1.1021 + if (_excess[u] < 0) {
1.1022 + _res_cap[a] = -_excess[u] + 1;
1.1023 + } else {
1.1024 + _res_cap[a] = 1;
1.1025 + }
1.1026 + _res_cap[_reverse[a]] = 0;
1.1027 + _cost[a] = 0;
1.1028 + _cost[_reverse[a]] = 0;
1.1029 + }
1.1030 + } else {
1.1031 + _pi[_root] = 0;
1.1032 + _excess[_root] = 0;
1.1033 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.1034 + _res_cap[a] = 1;
1.1035 + _res_cap[_reverse[a]] = 0;
1.1036 + _cost[a] = 0;
1.1037 + _cost[_reverse[a]] = 0;
1.1038 + }
1.1039 + }
1.1040 +
1.1041 + // Initialize delta value
1.1042 if (scaling) {
1.1043 // With scaling
1.1044 - Supply max_sup = 0, max_dem = 0;
1.1045 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.1046 - if ( _supply[n] > max_sup) max_sup = _supply[n];
1.1047 - if (-_supply[n] > max_dem) max_dem = -_supply[n];
1.1048 + Value max_sup = 0, max_dem = 0;
1.1049 + for (int i = 0; i != _node_num; ++i) {
1.1050 + if ( _excess[i] > max_sup) max_sup = _excess[i];
1.1051 + if (-_excess[i] > max_dem) max_dem = -_excess[i];
1.1052 }
1.1053 - Capacity max_cap = 0;
1.1054 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.1055 - if (_capacity[e] > max_cap) max_cap = _capacity[e];
1.1056 + Value max_cap = 0;
1.1057 + for (int j = 0; j != _res_arc_num; ++j) {
1.1058 + if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
1.1059 }
1.1060 max_sup = std::min(std::min(max_sup, max_dem), max_cap);
1.1061 _phase_num = 0;
1.1062 @@ -533,53 +716,70 @@
1.1063 _delta = 1;
1.1064 }
1.1065
1.1066 - return true;
1.1067 + return OPTIMAL;
1.1068 }
1.1069
1.1070 - bool start() {
1.1071 + ProblemType start() {
1.1072 + // Execute the algorithm
1.1073 + ProblemType pt;
1.1074 if (_delta > 1)
1.1075 - return startWithScaling();
1.1076 + pt = startWithScaling();
1.1077 else
1.1078 - return startWithoutScaling();
1.1079 + pt = startWithoutScaling();
1.1080 +
1.1081 + // Handle non-zero lower bounds
1.1082 + if (_have_lower) {
1.1083 + for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) {
1.1084 + if (!_forward[j]) _res_cap[j] += _lower[j];
1.1085 + }
1.1086 + }
1.1087 +
1.1088 + // Shift potentials if necessary
1.1089 + Cost pr = _pi[_root];
1.1090 + if (_sum_supply < 0 || pr > 0) {
1.1091 + for (int i = 0; i != _node_num; ++i) {
1.1092 + _pi[i] -= pr;
1.1093 + }
1.1094 + }
1.1095 +
1.1096 + return pt;
1.1097 }
1.1098
1.1099 - /// Execute the capacity scaling algorithm.
1.1100 - bool startWithScaling() {
1.1101 - // Processing capacity scaling phases
1.1102 - Node s, t;
1.1103 + // Execute the capacity scaling algorithm
1.1104 + ProblemType startWithScaling() {
1.1105 + // Process capacity scaling phases
1.1106 + int s, t;
1.1107 int phase_cnt = 0;
1.1108 int factor = 4;
1.1109 + ResidualDijkstra _dijkstra(*this);
1.1110 while (true) {
1.1111 - // Saturating all arcs not satisfying the optimality condition
1.1112 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.1113 - Node u = _graph.source(e), v = _graph.target(e);
1.1114 - Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
1.1115 - if (c < 0 && _res_cap[e] >= _delta) {
1.1116 - _excess[u] -= _res_cap[e];
1.1117 - _excess[v] += _res_cap[e];
1.1118 - (*_flow)[e] = _capacity[e];
1.1119 - _res_cap[e] = 0;
1.1120 - }
1.1121 - else if (c > 0 && (*_flow)[e] >= _delta) {
1.1122 - _excess[u] += (*_flow)[e];
1.1123 - _excess[v] -= (*_flow)[e];
1.1124 - (*_flow)[e] = 0;
1.1125 - _res_cap[e] = _capacity[e];
1.1126 + // Saturate all arcs not satisfying the optimality condition
1.1127 + for (int u = 0; u != _node_num; ++u) {
1.1128 + for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
1.1129 + int v = _target[a];
1.1130 + Cost c = _cost[a] + _pi[u] - _pi[v];
1.1131 + Value rc = _res_cap[a];
1.1132 + if (c < 0 && rc >= _delta) {
1.1133 + _excess[u] -= rc;
1.1134 + _excess[v] += rc;
1.1135 + _res_cap[a] = 0;
1.1136 + _res_cap[_reverse[a]] += rc;
1.1137 + }
1.1138 }
1.1139 }
1.1140
1.1141 - // Finding excess nodes and deficit nodes
1.1142 + // Find excess nodes and deficit nodes
1.1143 _excess_nodes.clear();
1.1144 _deficit_nodes.clear();
1.1145 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.1146 - if (_excess[n] >= _delta) _excess_nodes.push_back(n);
1.1147 - if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
1.1148 + for (int u = 0; u != _node_num; ++u) {
1.1149 + if (_excess[u] >= _delta) _excess_nodes.push_back(u);
1.1150 + if (_excess[u] <= -_delta) _deficit_nodes.push_back(u);
1.1151 }
1.1152 int next_node = 0, next_def_node = 0;
1.1153
1.1154 - // Finding augmenting shortest paths
1.1155 + // Find augmenting shortest paths
1.1156 while (next_node < int(_excess_nodes.size())) {
1.1157 - // Checking deficit nodes
1.1158 + // Check deficit nodes
1.1159 if (_delta > 1) {
1.1160 bool delta_deficit = false;
1.1161 for ( ; next_def_node < int(_deficit_nodes.size());
1.1162 @@ -592,44 +792,31 @@
1.1163 if (!delta_deficit) break;
1.1164 }
1.1165
1.1166 - // Running Dijkstra
1.1167 + // Run Dijkstra in the residual network
1.1168 s = _excess_nodes[next_node];
1.1169 - if ((t = _dijkstra->run(s, _delta)) == INVALID) {
1.1170 + if ((t = _dijkstra.run(s, _delta)) == -1) {
1.1171 if (_delta > 1) {
1.1172 ++next_node;
1.1173 continue;
1.1174 }
1.1175 - return false;
1.1176 + return INFEASIBLE;
1.1177 }
1.1178
1.1179 - // Augmenting along a shortest path from s to t.
1.1180 - Capacity d = std::min(_excess[s], -_excess[t]);
1.1181 - Node u = t;
1.1182 - Arc e;
1.1183 + // Augment along a shortest path from s to t
1.1184 + Value d = std::min(_excess[s], -_excess[t]);
1.1185 + int u = t;
1.1186 + int a;
1.1187 if (d > _delta) {
1.1188 - while ((e = _pred[u]) != INVALID) {
1.1189 - Capacity rc;
1.1190 - if (u == _graph.target(e)) {
1.1191 - rc = _res_cap[e];
1.1192 - u = _graph.source(e);
1.1193 - } else {
1.1194 - rc = (*_flow)[e];
1.1195 - u = _graph.target(e);
1.1196 - }
1.1197 - if (rc < d) d = rc;
1.1198 + while ((a = _pred[u]) != -1) {
1.1199 + if (_res_cap[a] < d) d = _res_cap[a];
1.1200 + u = _source[a];
1.1201 }
1.1202 }
1.1203 u = t;
1.1204 - while ((e = _pred[u]) != INVALID) {
1.1205 - if (u == _graph.target(e)) {
1.1206 - (*_flow)[e] += d;
1.1207 - _res_cap[e] -= d;
1.1208 - u = _graph.source(e);
1.1209 - } else {
1.1210 - (*_flow)[e] -= d;
1.1211 - _res_cap[e] += d;
1.1212 - u = _graph.target(e);
1.1213 - }
1.1214 + while ((a = _pred[u]) != -1) {
1.1215 + _res_cap[a] -= d;
1.1216 + _res_cap[_reverse[a]] += d;
1.1217 + u = _source[a];
1.1218 }
1.1219 _excess[s] -= d;
1.1220 _excess[t] += d;
1.1221 @@ -638,74 +825,54 @@
1.1222 }
1.1223
1.1224 if (_delta == 1) break;
1.1225 - if (++phase_cnt > _phase_num / 4) factor = 2;
1.1226 + if (++phase_cnt == _phase_num / 4) factor = 2;
1.1227 _delta = _delta <= factor ? 1 : _delta / factor;
1.1228 }
1.1229
1.1230 - // Handling non-zero lower bounds
1.1231 - if (_lower) {
1.1232 - for (ArcIt e(_graph); e != INVALID; ++e)
1.1233 - (*_flow)[e] += (*_lower)[e];
1.1234 - }
1.1235 - return true;
1.1236 + return OPTIMAL;
1.1237 }
1.1238
1.1239 - /// Execute the successive shortest path algorithm.
1.1240 - bool startWithoutScaling() {
1.1241 - // Finding excess nodes
1.1242 - for (NodeIt n(_graph); n != INVALID; ++n)
1.1243 - if (_excess[n] > 0) _excess_nodes.push_back(n);
1.1244 - if (_excess_nodes.size() == 0) return true;
1.1245 + // Execute the successive shortest path algorithm
1.1246 + ProblemType startWithoutScaling() {
1.1247 + // Find excess nodes
1.1248 + _excess_nodes.clear();
1.1249 + for (int i = 0; i != _node_num; ++i) {
1.1250 + if (_excess[i] > 0) _excess_nodes.push_back(i);
1.1251 + }
1.1252 + if (_excess_nodes.size() == 0) return OPTIMAL;
1.1253 int next_node = 0;
1.1254
1.1255 - // Finding shortest paths
1.1256 - Node s, t;
1.1257 + // Find shortest paths
1.1258 + int s, t;
1.1259 + ResidualDijkstra _dijkstra(*this);
1.1260 while ( _excess[_excess_nodes[next_node]] > 0 ||
1.1261 ++next_node < int(_excess_nodes.size()) )
1.1262 {
1.1263 - // Running Dijkstra
1.1264 + // Run Dijkstra in the residual network
1.1265 s = _excess_nodes[next_node];
1.1266 - if ((t = _dijkstra->run(s)) == INVALID) return false;
1.1267 + if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
1.1268
1.1269 - // Augmenting along a shortest path from s to t
1.1270 - Capacity d = std::min(_excess[s], -_excess[t]);
1.1271 - Node u = t;
1.1272 - Arc e;
1.1273 + // Augment along a shortest path from s to t
1.1274 + Value d = std::min(_excess[s], -_excess[t]);
1.1275 + int u = t;
1.1276 + int a;
1.1277 if (d > 1) {
1.1278 - while ((e = _pred[u]) != INVALID) {
1.1279 - Capacity rc;
1.1280 - if (u == _graph.target(e)) {
1.1281 - rc = _res_cap[e];
1.1282 - u = _graph.source(e);
1.1283 - } else {
1.1284 - rc = (*_flow)[e];
1.1285 - u = _graph.target(e);
1.1286 - }
1.1287 - if (rc < d) d = rc;
1.1288 + while ((a = _pred[u]) != -1) {
1.1289 + if (_res_cap[a] < d) d = _res_cap[a];
1.1290 + u = _source[a];
1.1291 }
1.1292 }
1.1293 u = t;
1.1294 - while ((e = _pred[u]) != INVALID) {
1.1295 - if (u == _graph.target(e)) {
1.1296 - (*_flow)[e] += d;
1.1297 - _res_cap[e] -= d;
1.1298 - u = _graph.source(e);
1.1299 - } else {
1.1300 - (*_flow)[e] -= d;
1.1301 - _res_cap[e] += d;
1.1302 - u = _graph.target(e);
1.1303 - }
1.1304 + while ((a = _pred[u]) != -1) {
1.1305 + _res_cap[a] -= d;
1.1306 + _res_cap[_reverse[a]] += d;
1.1307 + u = _source[a];
1.1308 }
1.1309 _excess[s] -= d;
1.1310 _excess[t] += d;
1.1311 }
1.1312
1.1313 - // Handling non-zero lower bounds
1.1314 - if (_lower) {
1.1315 - for (ArcIt e(_graph); e != INVALID; ++e)
1.1316 - (*_flow)[e] += (*_lower)[e];
1.1317 - }
1.1318 - return true;
1.1319 + return OPTIMAL;
1.1320 }
1.1321
1.1322 }; //class CapacityScaling