COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/capacity_scaling.h

1.2
Last change on this file was 1004:1e87c18cf65e, checked in by Alpar Juttner <alpar@…>, 4 years ago

Merge bugfix #600 to branch 1.2

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