COIN-OR::LEMON - Graph Library

source: lemon/lemon/capacity_scaling.h @ 872:fa6f37d7a25b

Last change on this file since 872:fa6f37d7a25b was 872:fa6f37d7a25b, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Entirely rework CapacityScaling? (#180)

  • Use the new interface similarly to NetworkSimplex?.
  • Rework the implementation using an efficient internal structure for handling the residual network. This improvement made the code much faster (up to 2-5 times faster on large graphs).
  • Handle GEQ supply type (LEQ is not supported).
  • Handle negative costs for arcs of finite capacity. (Note that this algorithm cannot handle arcs of negative cost and infinite upper bound, thus it returns UNBOUNDED if such an arc exists.)
  • Extend the documentation.
File size: 26.8 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
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.
12 *
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
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_CAPACITY_SCALING_H
20#define LEMON_CAPACITY_SCALING_H
21
22/// \ingroup min_cost_flow_algs
23///
24/// \file
25/// \brief Capacity Scaling algorithm for finding a minimum cost flow.
26
27#include <vector>
28#include <limits>
29#include <lemon/core.h>
30#include <lemon/bin_heap.h>
31
32namespace lemon {
33
34  /// \addtogroup min_cost_flow_algs
35  /// @{
36
37  /// \brief Implementation of the Capacity Scaling algorithm for
38  /// finding a \ref min_cost_flow "minimum cost flow".
39  ///
40  /// \ref CapacityScaling implements the capacity scaling version
41  /// of the successive shortest path algorithm for finding a
42  /// \ref min_cost_flow "minimum cost flow". It is an efficient dual
43  /// solution method.
44  ///
45  /// Most of the parameters of the problem (except for the digraph)
46  /// can be given using separate functions, and the algorithm can be
47  /// executed using the \ref run() function. If some parameters are not
48  /// specified, then default values will be used.
49  ///
50  /// \tparam GR The digraph type the algorithm runs on.
51  /// \tparam V The value type used for flow amounts, capacity bounds
52  /// and supply values in the algorithm. By default it is \c int.
53  /// \tparam C The value type used for costs and potentials in the
54  /// algorithm. By default it is the same as \c V.
55  ///
56  /// \warning Both value types must be signed and all input data must
57  /// be integer.
58  /// \warning This algorithm does not support negative costs for such
59  /// arcs that have infinite upper bound.
60  template <typename GR, typename V = int, typename C = V>
61  class CapacityScaling
62  {
63  public:
64
65    /// The type of the flow amounts, capacity bounds and supply values
66    typedef V Value;
67    /// The type of the arc costs
68    typedef C Cost;
69
70  public:
71
72    /// \brief Problem type constants for the \c run() function.
73    ///
74    /// Enum type containing the problem type constants that can be
75    /// returned by the \ref run() function of the algorithm.
76    enum ProblemType {
77      /// The problem has no feasible solution (flow).
78      INFEASIBLE,
79      /// The problem has optimal solution (i.e. it is feasible and
80      /// bounded), and the algorithm has found optimal flow and node
81      /// potentials (primal and dual solutions).
82      OPTIMAL,
83      /// The digraph contains an arc of negative cost and infinite
84      /// upper bound. It means that the objective function is unbounded
85      /// on that arc, however note that it could actually be bounded
86      /// over the feasible flows, but this algroithm cannot handle
87      /// these cases.
88      UNBOUNDED
89    };
90 
91  private:
92
93    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
94
95    typedef std::vector<Arc> ArcVector;
96    typedef std::vector<Node> NodeVector;
97    typedef std::vector<int> IntVector;
98    typedef std::vector<bool> BoolVector;
99    typedef std::vector<Value> ValueVector;
100    typedef std::vector<Cost> CostVector;
101
102  private:
103
104    // Data related to the underlying digraph
105    const GR &_graph;
106    int _node_num;
107    int _arc_num;
108    int _res_arc_num;
109    int _root;
110
111    // Parameters of the problem
112    bool _have_lower;
113    Value _sum_supply;
114
115    // Data structures for storing the digraph
116    IntNodeMap _node_id;
117    IntArcMap _arc_idf;
118    IntArcMap _arc_idb;
119    IntVector _first_out;
120    BoolVector _forward;
121    IntVector _source;
122    IntVector _target;
123    IntVector _reverse;
124
125    // Node and arc data
126    ValueVector _lower;
127    ValueVector _upper;
128    CostVector _cost;
129    ValueVector _supply;
130
131    ValueVector _res_cap;
132    CostVector _pi;
133    ValueVector _excess;
134    IntVector _excess_nodes;
135    IntVector _deficit_nodes;
136
137    Value _delta;
138    int _phase_num;
139    IntVector _pred;
140
141  public:
142 
143    /// \brief Constant for infinite upper bounds (capacities).
144    ///
145    /// Constant for infinite upper bounds (capacities).
146    /// It is \c std::numeric_limits<Value>::infinity() if available,
147    /// \c std::numeric_limits<Value>::max() otherwise.
148    const Value INF;
149
150  private:
151
152    // Special implementation of the Dijkstra algorithm for finding
153    // shortest paths in the residual network of the digraph with
154    // respect to the reduced arc costs and modifying the node
155    // potentials according to the found distance labels.
156    class ResidualDijkstra
157    {
158      typedef RangeMap<int> HeapCrossRef;
159      typedef BinHeap<Cost, HeapCrossRef> Heap;
160
161    private:
162
163      int _node_num;
164      const IntVector &_first_out;
165      const IntVector &_target;
166      const CostVector &_cost;
167      const ValueVector &_res_cap;
168      const ValueVector &_excess;
169      CostVector &_pi;
170      IntVector &_pred;
171     
172      IntVector _proc_nodes;
173      CostVector _dist;
174     
175    public:
176
177      ResidualDijkstra(CapacityScaling& cs) :
178        _node_num(cs._node_num), _first_out(cs._first_out),
179        _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap),
180        _excess(cs._excess), _pi(cs._pi), _pred(cs._pred),
181        _dist(cs._node_num)
182      {}
183
184      int run(int s, Value delta = 1) {
185        HeapCrossRef heap_cross_ref(_node_num, Heap::PRE_HEAP);
186        Heap heap(heap_cross_ref);
187        heap.push(s, 0);
188        _pred[s] = -1;
189        _proc_nodes.clear();
190
191        // Process nodes
192        while (!heap.empty() && _excess[heap.top()] > -delta) {
193          int u = heap.top(), v;
194          Cost d = heap.prio() + _pi[u], dn;
195          _dist[u] = heap.prio();
196          _proc_nodes.push_back(u);
197          heap.pop();
198
199          // Traverse outgoing residual arcs
200          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
201            if (_res_cap[a] < delta) continue;
202            v = _target[a];
203            switch (heap.state(v)) {
204              case Heap::PRE_HEAP:
205                heap.push(v, d + _cost[a] - _pi[v]);
206                _pred[v] = a;
207                break;
208              case Heap::IN_HEAP:
209                dn = d + _cost[a] - _pi[v];
210                if (dn < heap[v]) {
211                  heap.decrease(v, dn);
212                  _pred[v] = a;
213                }
214                break;
215              case Heap::POST_HEAP:
216                break;
217            }
218          }
219        }
220        if (heap.empty()) return -1;
221
222        // Update potentials of processed nodes
223        int t = heap.top();
224        Cost dt = heap.prio();
225        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
226          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
227        }
228
229        return t;
230      }
231
232    }; //class ResidualDijkstra
233
234  public:
235
236    /// \brief Constructor.
237    ///
238    /// The constructor of the class.
239    ///
240    /// \param graph The digraph the algorithm runs on.
241    CapacityScaling(const GR& graph) :
242      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
243      INF(std::numeric_limits<Value>::has_infinity ?
244          std::numeric_limits<Value>::infinity() :
245          std::numeric_limits<Value>::max())
246    {
247      // Check the value types
248      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
249        "The flow type of CapacityScaling must be signed");
250      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
251        "The cost type of CapacityScaling must be signed");
252
253      // Resize vectors
254      _node_num = countNodes(_graph);
255      _arc_num = countArcs(_graph);
256      _res_arc_num = 2 * (_arc_num + _node_num);
257      _root = _node_num;
258      ++_node_num;
259
260      _first_out.resize(_node_num + 1);
261      _forward.resize(_res_arc_num);
262      _source.resize(_res_arc_num);
263      _target.resize(_res_arc_num);
264      _reverse.resize(_res_arc_num);
265
266      _lower.resize(_res_arc_num);
267      _upper.resize(_res_arc_num);
268      _cost.resize(_res_arc_num);
269      _supply.resize(_node_num);
270     
271      _res_cap.resize(_res_arc_num);
272      _pi.resize(_node_num);
273      _excess.resize(_node_num);
274      _pred.resize(_node_num);
275
276      // Copy the graph
277      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
278      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
279        _node_id[n] = i;
280      }
281      i = 0;
282      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
283        _first_out[i] = j;
284        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
285          _arc_idf[a] = j;
286          _forward[j] = true;
287          _source[j] = i;
288          _target[j] = _node_id[_graph.runningNode(a)];
289        }
290        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
291          _arc_idb[a] = j;
292          _forward[j] = false;
293          _source[j] = i;
294          _target[j] = _node_id[_graph.runningNode(a)];
295        }
296        _forward[j] = false;
297        _source[j] = i;
298        _target[j] = _root;
299        _reverse[j] = k;
300        _forward[k] = true;
301        _source[k] = _root;
302        _target[k] = i;
303        _reverse[k] = j;
304        ++j; ++k;
305      }
306      _first_out[i] = j;
307      _first_out[_node_num] = k;
308      for (ArcIt a(_graph); a != INVALID; ++a) {
309        int fi = _arc_idf[a];
310        int bi = _arc_idb[a];
311        _reverse[fi] = bi;
312        _reverse[bi] = fi;
313      }
314     
315      // Reset parameters
316      reset();
317    }
318
319    /// \name Parameters
320    /// The parameters of the algorithm can be specified using these
321    /// functions.
322
323    /// @{
324
325    /// \brief Set the lower bounds on the arcs.
326    ///
327    /// This function sets the lower bounds on the arcs.
328    /// If it is not used before calling \ref run(), the lower bounds
329    /// will be set to zero on all arcs.
330    ///
331    /// \param map An arc map storing the lower bounds.
332    /// Its \c Value type must be convertible to the \c Value type
333    /// of the algorithm.
334    ///
335    /// \return <tt>(*this)</tt>
336    template <typename LowerMap>
337    CapacityScaling& lowerMap(const LowerMap& map) {
338      _have_lower = true;
339      for (ArcIt a(_graph); a != INVALID; ++a) {
340        _lower[_arc_idf[a]] = map[a];
341        _lower[_arc_idb[a]] = map[a];
342      }
343      return *this;
344    }
345
346    /// \brief Set the upper bounds (capacities) on the arcs.
347    ///
348    /// This function sets the upper bounds (capacities) on the arcs.
349    /// If it is not used before calling \ref run(), the upper bounds
350    /// will be set to \ref INF on all arcs (i.e. the flow value will be
351    /// unbounded from above on each arc).
352    ///
353    /// \param map An arc map storing the upper bounds.
354    /// Its \c Value type must be convertible to the \c Value type
355    /// of the algorithm.
356    ///
357    /// \return <tt>(*this)</tt>
358    template<typename UpperMap>
359    CapacityScaling& upperMap(const UpperMap& map) {
360      for (ArcIt a(_graph); a != INVALID; ++a) {
361        _upper[_arc_idf[a]] = map[a];
362      }
363      return *this;
364    }
365
366    /// \brief Set the costs of the arcs.
367    ///
368    /// This function sets the costs of the arcs.
369    /// If it is not used before calling \ref run(), the costs
370    /// will be set to \c 1 on all arcs.
371    ///
372    /// \param map An arc map storing the costs.
373    /// Its \c Value type must be convertible to the \c Cost type
374    /// of the algorithm.
375    ///
376    /// \return <tt>(*this)</tt>
377    template<typename CostMap>
378    CapacityScaling& costMap(const CostMap& map) {
379      for (ArcIt a(_graph); a != INVALID; ++a) {
380        _cost[_arc_idf[a]] =  map[a];
381        _cost[_arc_idb[a]] = -map[a];
382      }
383      return *this;
384    }
385
386    /// \brief Set the supply values of the nodes.
387    ///
388    /// This function sets the supply values of the nodes.
389    /// If neither this function nor \ref stSupply() is used before
390    /// calling \ref run(), the supply of each node will be set to zero.
391    ///
392    /// \param map A node map storing the supply values.
393    /// Its \c Value type must be convertible to the \c Value type
394    /// of the algorithm.
395    ///
396    /// \return <tt>(*this)</tt>
397    template<typename SupplyMap>
398    CapacityScaling& supplyMap(const SupplyMap& map) {
399      for (NodeIt n(_graph); n != INVALID; ++n) {
400        _supply[_node_id[n]] = map[n];
401      }
402      return *this;
403    }
404
405    /// \brief Set single source and target nodes and a supply value.
406    ///
407    /// This function sets a single source node and a single target node
408    /// and the required flow value.
409    /// If neither this function nor \ref supplyMap() is used before
410    /// calling \ref run(), the supply of each node will be set to zero.
411    ///
412    /// Using this function has the same effect as using \ref supplyMap()
413    /// with such a map in which \c k is assigned to \c s, \c -k is
414    /// assigned to \c t and all other nodes have zero supply value.
415    ///
416    /// \param s The source node.
417    /// \param t The target node.
418    /// \param k The required amount of flow from node \c s to node \c t
419    /// (i.e. the supply of \c s and the demand of \c t).
420    ///
421    /// \return <tt>(*this)</tt>
422    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
423      for (int i = 0; i != _node_num; ++i) {
424        _supply[i] = 0;
425      }
426      _supply[_node_id[s]] =  k;
427      _supply[_node_id[t]] = -k;
428      return *this;
429    }
430   
431    /// @}
432
433    /// \name Execution control
434
435    /// @{
436
437    /// \brief Run the algorithm.
438    ///
439    /// This function runs the algorithm.
440    /// The paramters can be specified using functions \ref lowerMap(),
441    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
442    /// For example,
443    /// \code
444    ///   CapacityScaling<ListDigraph> cs(graph);
445    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
446    ///     .supplyMap(sup).run();
447    /// \endcode
448    ///
449    /// This function can be called more than once. All the parameters
450    /// that have been given are kept for the next call, unless
451    /// \ref reset() is called, thus only the modified parameters
452    /// have to be set again. See \ref reset() for examples.
453    /// However the underlying digraph must not be modified after this
454    /// class have been constructed, since it copies the digraph.
455    ///
456    /// \param scaling Enable or disable capacity scaling.
457    /// If the maximum upper bound and/or the amount of total supply
458    /// is rather small, the algorithm could be slightly faster without
459    /// scaling.
460    ///
461    /// \return \c INFEASIBLE if no feasible flow exists,
462    /// \n \c OPTIMAL if the problem has optimal solution
463    /// (i.e. it is feasible and bounded), and the algorithm has found
464    /// optimal flow and node potentials (primal and dual solutions),
465    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
466    /// and infinite upper bound. It means that the objective function
467    /// is unbounded on that arc, however note that it could actually be
468    /// bounded over the feasible flows, but this algroithm cannot handle
469    /// these cases.
470    ///
471    /// \see ProblemType
472    ProblemType run(bool scaling = true) {
473      ProblemType pt = init(scaling);
474      if (pt != OPTIMAL) return pt;
475      return start();
476    }
477
478    /// \brief Reset all the parameters that have been given before.
479    ///
480    /// This function resets all the paramaters that have been given
481    /// before using functions \ref lowerMap(), \ref upperMap(),
482    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
483    ///
484    /// It is useful for multiple run() calls. If this function is not
485    /// used, all the parameters given before are kept for the next
486    /// \ref run() call.
487    /// However the underlying digraph must not be modified after this
488    /// class have been constructed, since it copies and extends the graph.
489    ///
490    /// For example,
491    /// \code
492    ///   CapacityScaling<ListDigraph> cs(graph);
493    ///
494    ///   // First run
495    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
496    ///     .supplyMap(sup).run();
497    ///
498    ///   // Run again with modified cost map (reset() is not called,
499    ///   // so only the cost map have to be set again)
500    ///   cost[e] += 100;
501    ///   cs.costMap(cost).run();
502    ///
503    ///   // Run again from scratch using reset()
504    ///   // (the lower bounds will be set to zero on all arcs)
505    ///   cs.reset();
506    ///   cs.upperMap(capacity).costMap(cost)
507    ///     .supplyMap(sup).run();
508    /// \endcode
509    ///
510    /// \return <tt>(*this)</tt>
511    CapacityScaling& reset() {
512      for (int i = 0; i != _node_num; ++i) {
513        _supply[i] = 0;
514      }
515      for (int j = 0; j != _res_arc_num; ++j) {
516        _lower[j] = 0;
517        _upper[j] = INF;
518        _cost[j] = _forward[j] ? 1 : -1;
519      }
520      _have_lower = false;
521      return *this;
522    }
523
524    /// @}
525
526    /// \name Query Functions
527    /// The results of the algorithm can be obtained using these
528    /// functions.\n
529    /// The \ref run() function must be called before using them.
530
531    /// @{
532
533    /// \brief Return the total cost of the found flow.
534    ///
535    /// This function returns the total cost of the found flow.
536    /// Its complexity is O(e).
537    ///
538    /// \note The return type of the function can be specified as a
539    /// template parameter. For example,
540    /// \code
541    ///   cs.totalCost<double>();
542    /// \endcode
543    /// It is useful if the total cost cannot be stored in the \c Cost
544    /// type of the algorithm, which is the default return type of the
545    /// function.
546    ///
547    /// \pre \ref run() must be called before using this function.
548    template <typename Number>
549    Number totalCost() const {
550      Number c = 0;
551      for (ArcIt a(_graph); a != INVALID; ++a) {
552        int i = _arc_idb[a];
553        c += static_cast<Number>(_res_cap[i]) *
554             (-static_cast<Number>(_cost[i]));
555      }
556      return c;
557    }
558
559#ifndef DOXYGEN
560    Cost totalCost() const {
561      return totalCost<Cost>();
562    }
563#endif
564
565    /// \brief Return the flow on the given arc.
566    ///
567    /// This function returns the flow on the given arc.
568    ///
569    /// \pre \ref run() must be called before using this function.
570    Value flow(const Arc& a) const {
571      return _res_cap[_arc_idb[a]];
572    }
573
574    /// \brief Return the flow map (the primal solution).
575    ///
576    /// This function copies the flow value on each arc into the given
577    /// map. The \c Value type of the algorithm must be convertible to
578    /// the \c Value type of the map.
579    ///
580    /// \pre \ref run() must be called before using this function.
581    template <typename FlowMap>
582    void flowMap(FlowMap &map) const {
583      for (ArcIt a(_graph); a != INVALID; ++a) {
584        map.set(a, _res_cap[_arc_idb[a]]);
585      }
586    }
587
588    /// \brief Return the potential (dual value) of the given node.
589    ///
590    /// This function returns the potential (dual value) of the
591    /// given node.
592    ///
593    /// \pre \ref run() must be called before using this function.
594    Cost potential(const Node& n) const {
595      return _pi[_node_id[n]];
596    }
597
598    /// \brief Return the potential map (the dual solution).
599    ///
600    /// This function copies the potential (dual value) of each node
601    /// into the given map.
602    /// The \c Cost type of the algorithm must be convertible to the
603    /// \c Value type of the map.
604    ///
605    /// \pre \ref run() must be called before using this function.
606    template <typename PotentialMap>
607    void potentialMap(PotentialMap &map) const {
608      for (NodeIt n(_graph); n != INVALID; ++n) {
609        map.set(n, _pi[_node_id[n]]);
610      }
611    }
612
613    /// @}
614
615  private:
616
617    // Initialize the algorithm
618    ProblemType init(bool scaling) {
619      if (_node_num == 0) return INFEASIBLE;
620
621      // Check the sum of supply values
622      _sum_supply = 0;
623      for (int i = 0; i != _root; ++i) {
624        _sum_supply += _supply[i];
625      }
626      if (_sum_supply > 0) return INFEASIBLE;
627     
628      // Initialize maps
629      for (int i = 0; i != _root; ++i) {
630        _pi[i] = 0;
631        _excess[i] = _supply[i];
632      }
633
634      // Remove non-zero lower bounds
635      if (_have_lower) {
636        for (int i = 0; i != _root; ++i) {
637          for (int j = _first_out[i]; j != _first_out[i+1]; ++j) {
638            if (_forward[j]) {
639              Value c = _lower[j];
640              if (c >= 0) {
641                _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF;
642              } else {
643                _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF;
644              }
645              _excess[i] -= c;
646              _excess[_target[j]] += c;
647            } else {
648              _res_cap[j] = 0;
649            }
650          }
651        }
652      } else {
653        for (int j = 0; j != _res_arc_num; ++j) {
654          _res_cap[j] = _forward[j] ? _upper[j] : 0;
655        }
656      }
657
658      // Handle negative costs
659      for (int u = 0; u != _root; ++u) {
660        for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
661          Value rc = _res_cap[a];
662          if (_cost[a] < 0 && rc > 0) {
663            if (rc == INF) return UNBOUNDED;
664            _excess[u] -= rc;
665            _excess[_target[a]] += rc;
666            _res_cap[a] = 0;
667            _res_cap[_reverse[a]] += rc;
668          }
669        }
670      }
671     
672      // Handle GEQ supply type
673      if (_sum_supply < 0) {
674        _pi[_root] = 0;
675        _excess[_root] = -_sum_supply;
676        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
677          int u = _target[a];
678          if (_excess[u] < 0) {
679            _res_cap[a] = -_excess[u] + 1;
680          } else {
681            _res_cap[a] = 1;
682          }
683          _res_cap[_reverse[a]] = 0;
684          _cost[a] = 0;
685          _cost[_reverse[a]] = 0;
686        }
687      } else {
688        _pi[_root] = 0;
689        _excess[_root] = 0;
690        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
691          _res_cap[a] = 1;
692          _res_cap[_reverse[a]] = 0;
693          _cost[a] = 0;
694          _cost[_reverse[a]] = 0;
695        }
696      }
697
698      // Initialize delta value
699      if (scaling) {
700        // With scaling
701        Value max_sup = 0, max_dem = 0;
702        for (int i = 0; i != _node_num; ++i) {
703          if ( _excess[i] > max_sup) max_sup =  _excess[i];
704          if (-_excess[i] > max_dem) max_dem = -_excess[i];
705        }
706        Value max_cap = 0;
707        for (int j = 0; j != _res_arc_num; ++j) {
708          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
709        }
710        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
711        _phase_num = 0;
712        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
713          ++_phase_num;
714      } else {
715        // Without scaling
716        _delta = 1;
717      }
718
719      return OPTIMAL;
720    }
721
722    ProblemType start() {
723      // Execute the algorithm
724      ProblemType pt;
725      if (_delta > 1)
726        pt = startWithScaling();
727      else
728        pt = startWithoutScaling();
729
730      // Handle non-zero lower bounds
731      if (_have_lower) {
732        for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) {
733          if (!_forward[j]) _res_cap[j] += _lower[j];
734        }
735      }
736
737      // Shift potentials if necessary
738      Cost pr = _pi[_root];
739      if (_sum_supply < 0 || pr > 0) {
740        for (int i = 0; i != _node_num; ++i) {
741          _pi[i] -= pr;
742        }       
743      }
744     
745      return pt;
746    }
747
748    // Execute the capacity scaling algorithm
749    ProblemType startWithScaling() {
750      // Process capacity scaling phases
751      int s, t;
752      int phase_cnt = 0;
753      int factor = 4;
754      ResidualDijkstra _dijkstra(*this);
755      while (true) {
756        // Saturate all arcs not satisfying the optimality condition
757        for (int u = 0; u != _node_num; ++u) {
758          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
759            int v = _target[a];
760            Cost c = _cost[a] + _pi[u] - _pi[v];
761            Value rc = _res_cap[a];
762            if (c < 0 && rc >= _delta) {
763              _excess[u] -= rc;
764              _excess[v] += rc;
765              _res_cap[a] = 0;
766              _res_cap[_reverse[a]] += rc;
767            }
768          }
769        }
770
771        // Find excess nodes and deficit nodes
772        _excess_nodes.clear();
773        _deficit_nodes.clear();
774        for (int u = 0; u != _node_num; ++u) {
775          if (_excess[u] >=  _delta) _excess_nodes.push_back(u);
776          if (_excess[u] <= -_delta) _deficit_nodes.push_back(u);
777        }
778        int next_node = 0, next_def_node = 0;
779
780        // Find augmenting shortest paths
781        while (next_node < int(_excess_nodes.size())) {
782          // Check deficit nodes
783          if (_delta > 1) {
784            bool delta_deficit = false;
785            for ( ; next_def_node < int(_deficit_nodes.size());
786                    ++next_def_node ) {
787              if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
788                delta_deficit = true;
789                break;
790              }
791            }
792            if (!delta_deficit) break;
793          }
794
795          // Run Dijkstra in the residual network
796          s = _excess_nodes[next_node];
797          if ((t = _dijkstra.run(s, _delta)) == -1) {
798            if (_delta > 1) {
799              ++next_node;
800              continue;
801            }
802            return INFEASIBLE;
803          }
804
805          // Augment along a shortest path from s to t
806          Value d = std::min(_excess[s], -_excess[t]);
807          int u = t;
808          int a;
809          if (d > _delta) {
810            while ((a = _pred[u]) != -1) {
811              if (_res_cap[a] < d) d = _res_cap[a];
812              u = _source[a];
813            }
814          }
815          u = t;
816          while ((a = _pred[u]) != -1) {
817            _res_cap[a] -= d;
818            _res_cap[_reverse[a]] += d;
819            u = _source[a];
820          }
821          _excess[s] -= d;
822          _excess[t] += d;
823
824          if (_excess[s] < _delta) ++next_node;
825        }
826
827        if (_delta == 1) break;
828        if (++phase_cnt == _phase_num / 4) factor = 2;
829        _delta = _delta <= factor ? 1 : _delta / factor;
830      }
831
832      return OPTIMAL;
833    }
834
835    // Execute the successive shortest path algorithm
836    ProblemType startWithoutScaling() {
837      // Find excess nodes
838      _excess_nodes.clear();
839      for (int i = 0; i != _node_num; ++i) {
840        if (_excess[i] > 0) _excess_nodes.push_back(i);
841      }
842      if (_excess_nodes.size() == 0) return OPTIMAL;
843      int next_node = 0;
844
845      // Find shortest paths
846      int s, t;
847      ResidualDijkstra _dijkstra(*this);
848      while ( _excess[_excess_nodes[next_node]] > 0 ||
849              ++next_node < int(_excess_nodes.size()) )
850      {
851        // Run Dijkstra in the residual network
852        s = _excess_nodes[next_node];
853        if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
854
855        // Augment along a shortest path from s to t
856        Value d = std::min(_excess[s], -_excess[t]);
857        int u = t;
858        int a;
859        if (d > 1) {
860          while ((a = _pred[u]) != -1) {
861            if (_res_cap[a] < d) d = _res_cap[a];
862            u = _source[a];
863          }
864        }
865        u = t;
866        while ((a = _pred[u]) != -1) {
867          _res_cap[a] -= d;
868          _res_cap[_reverse[a]] += d;
869          u = _source[a];
870        }
871        _excess[s] -= d;
872        _excess[t] += d;
873      }
874
875      return OPTIMAL;
876    }
877
878  }; //class CapacityScaling
879
880  ///@}
881
882} //namespace lemon
883
884#endif //LEMON_CAPACITY_SCALING_H
Note: See TracBrowser for help on using the repository browser.