1.1 --- a/lemon/cost_scaling.h	Tue Mar 15 19:32:21 2011 +0100
     1.2 +++ b/lemon/cost_scaling.h	Tue Mar 15 19:52:31 2011 +0100
     1.3 @@ -237,7 +237,6 @@
     1.4        std::vector<Value>& _v;
     1.5      };
     1.6  
     1.7 -    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
     1.8      typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
     1.9  
    1.10    private:
    1.11 @@ -288,14 +287,6 @@
    1.12      IntVector _rank;
    1.13      int _max_rank;
    1.14  
    1.15 -    // Data for a StaticDigraph structure
    1.16 -    typedef std::pair<int, int> IntPair;
    1.17 -    StaticDigraph _sgr;
    1.18 -    std::vector<IntPair> _arc_vec;
    1.19 -    std::vector<LargeCost> _cost_vec;
    1.20 -    LargeCostArcMap _cost_map;
    1.21 -    LargeCostNodeMap _pi_map;
    1.22 -
    1.23    public:
    1.24  
    1.25      /// \brief Constant for infinite upper bounds (capacities).
    1.26 @@ -342,7 +333,6 @@
    1.27      /// \param graph The digraph the algorithm runs on.
    1.28      CostScaling(const GR& graph) :
    1.29        _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
    1.30 -      _cost_map(_cost_vec), _pi_map(_pi),
    1.31        INF(std::numeric_limits<Value>::has_infinity ?
    1.32            std::numeric_limits<Value>::infinity() :
    1.33            std::numeric_limits<Value>::max())
    1.34 @@ -619,9 +609,6 @@
    1.35        _excess.resize(_res_node_num);
    1.36        _next_out.resize(_res_node_num);
    1.37  
    1.38 -      _arc_vec.reserve(_res_arc_num);
    1.39 -      _cost_vec.reserve(_res_arc_num);
    1.40 -
    1.41        // Copy the graph
    1.42        int i = 0, j = 0, k = 2 * _arc_num + _node_num;
    1.43        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    1.44 @@ -923,22 +910,63 @@
    1.45            break;
    1.46        }
    1.47  
    1.48 -      // Compute node potentials for the original costs
    1.49 -      _arc_vec.clear();
    1.50 -      _cost_vec.clear();
    1.51 -      for (int j = 0; j != _res_arc_num; ++j) {
    1.52 -        if (_res_cap[j] > 0) {
    1.53 -          _arc_vec.push_back(IntPair(_source[j], _target[j]));
    1.54 -          _cost_vec.push_back(_scost[j]);
    1.55 +      // Compute node potentials (dual solution)
    1.56 +      for (int i = 0; i != _res_node_num; ++i) {
    1.57 +        _pi[i] = static_cast<Cost>(_pi[i] / (_res_node_num * _alpha));
    1.58 +      }
    1.59 +      bool optimal = true;
    1.60 +      for (int i = 0; optimal && i != _res_node_num; ++i) {
    1.61 +        LargeCost pi_i = _pi[i];
    1.62 +        int last_out = _first_out[i+1];
    1.63 +        for (int j = _first_out[i]; j != last_out; ++j) {
    1.64 +          if (_res_cap[j] > 0 && _scost[j] + pi_i - _pi[_target[j]] < 0) {
    1.65 +            optimal = false;
    1.66 +            break;
    1.67 +          }
    1.68          }
    1.69        }
    1.70 -      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    1.71  
    1.72 -      typename BellmanFord<StaticDigraph, LargeCostArcMap>
    1.73 -        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
    1.74 -      bf.distMap(_pi_map);
    1.75 -      bf.init(0);
    1.76 -      bf.start();
    1.77 +      if (!optimal) {
    1.78 +        // Compute node potentials for the original costs with BellmanFord
    1.79 +        // (if it is necessary)
    1.80 +        typedef std::pair<int, int> IntPair;
    1.81 +        StaticDigraph sgr;
    1.82 +        std::vector<IntPair> arc_vec;
    1.83 +        std::vector<LargeCost> cost_vec;
    1.84 +        LargeCostArcMap cost_map(cost_vec);
    1.85 +
    1.86 +        arc_vec.clear();
    1.87 +        cost_vec.clear();
    1.88 +        for (int j = 0; j != _res_arc_num; ++j) {
    1.89 +          if (_res_cap[j] > 0) {
    1.90 +            int u = _source[j], v = _target[j];
    1.91 +            arc_vec.push_back(IntPair(u, v));
    1.92 +            cost_vec.push_back(_scost[j] + _pi[u] - _pi[v]);
    1.93 +          }
    1.94 +        }
    1.95 +        sgr.build(_res_node_num, arc_vec.begin(), arc_vec.end());
    1.96 +
    1.97 +        typename BellmanFord<StaticDigraph, LargeCostArcMap>::Create
    1.98 +          bf(sgr, cost_map);
    1.99 +        bf.init(0);
   1.100 +        bf.start();
   1.101 +
   1.102 +        for (int i = 0; i != _res_node_num; ++i) {
   1.103 +          _pi[i] += bf.dist(sgr.node(i));
   1.104 +        }
   1.105 +      }
   1.106 +
   1.107 +      // Shift potentials to meet the requirements of the GEQ type
   1.108 +      // optimality conditions
   1.109 +      LargeCost max_pot = _pi[_root];
   1.110 +      for (int i = 0; i != _res_node_num; ++i) {
   1.111 +        if (_pi[i] > max_pot) max_pot = _pi[i];
   1.112 +      }
   1.113 +      if (max_pot != 0) {
   1.114 +        for (int i = 0; i != _res_node_num; ++i) {
   1.115 +          _pi[i] -= max_pot;
   1.116 +        }
   1.117 +      }
   1.118  
   1.119        // Handle non-zero lower bounds
   1.120        if (_have_lower) {