Less map copying in NetworkSimplex (#234)
authorPeter Kovacs <kpeter@inf.elte.hu>
Wed, 29 Apr 2009 16:54:27 +0200
changeset 642111698359429
parent 641 756a5ec551c8
child 643 f3792d5bb294
Less map copying in NetworkSimplex (#234)

- The graph is copied in the constructor instead of the init() function.
It must not be modified after the class is constructed.
- The maps are copied once (instead of twice).
- Remove FlowMap, PotentialMap typedefs and flowMap(), pontentialMap()
setter functions.
- flowMap() and potentialMap() query functions copy the values into the
given map (reference) instead of returning a const reference to a
previously constructed map.
lemon/network_simplex.h
test/min_cost_flow_test.cc
     1.1 --- a/lemon/network_simplex.h	Wed Apr 29 14:25:51 2009 +0200
     1.2 +++ b/lemon/network_simplex.h	Wed Apr 29 16:54:27 2009 +0200
     1.3 @@ -72,21 +72,10 @@
     1.4    {
     1.5    public:
     1.6  
     1.7 -    /// The flow type of the algorithm
     1.8 +    /// The type of the flow amounts, capacity bounds and supply values
     1.9      typedef V Value;
    1.10 -    /// The cost type of the algorithm
    1.11 +    /// The type of the arc costs
    1.12      typedef C Cost;
    1.13 -#ifdef DOXYGEN
    1.14 -    /// The type of the flow map
    1.15 -    typedef GR::ArcMap<Value> FlowMap;
    1.16 -    /// The type of the potential map
    1.17 -    typedef GR::NodeMap<Cost> PotentialMap;
    1.18 -#else
    1.19 -    /// The type of the flow map
    1.20 -    typedef typename GR::template ArcMap<Value> FlowMap;
    1.21 -    /// The type of the potential map
    1.22 -    typedef typename GR::template NodeMap<Cost> PotentialMap;
    1.23 -#endif
    1.24  
    1.25    public:
    1.26  
    1.27 @@ -206,15 +195,11 @@
    1.28  
    1.29      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    1.30  
    1.31 -    typedef typename GR::template ArcMap<Value> ValueArcMap;
    1.32 -    typedef typename GR::template ArcMap<Cost> CostArcMap;
    1.33 -    typedef typename GR::template NodeMap<Value> ValueNodeMap;
    1.34 -
    1.35      typedef std::vector<Arc> ArcVector;
    1.36      typedef std::vector<Node> NodeVector;
    1.37      typedef std::vector<int> IntVector;
    1.38      typedef std::vector<bool> BoolVector;
    1.39 -    typedef std::vector<Value> FlowVector;
    1.40 +    typedef std::vector<Value> ValueVector;
    1.41      typedef std::vector<Cost> CostVector;
    1.42  
    1.43      // State constants for arcs
    1.44 @@ -232,34 +217,23 @@
    1.45      int _arc_num;
    1.46  
    1.47      // Parameters of the problem
    1.48 -    ValueArcMap *_plower;
    1.49 -    ValueArcMap *_pupper;
    1.50 -    CostArcMap *_pcost;
    1.51 -    ValueNodeMap *_psupply;
    1.52 -    bool _pstsup;
    1.53 -    Node _psource, _ptarget;
    1.54 -    Value _pstflow;
    1.55 +    bool _have_lower;
    1.56      SupplyType _stype;
    1.57 -    
    1.58      Value _sum_supply;
    1.59  
    1.60 -    // Result maps
    1.61 -    FlowMap *_flow_map;
    1.62 -    PotentialMap *_potential_map;
    1.63 -    bool _local_flow;
    1.64 -    bool _local_potential;
    1.65 -
    1.66      // Data structures for storing the digraph
    1.67      IntNodeMap _node_id;
    1.68 -    ArcVector _arc_ref;
    1.69 +    IntArcMap _arc_id;
    1.70      IntVector _source;
    1.71      IntVector _target;
    1.72  
    1.73      // Node and arc data
    1.74 -    FlowVector _cap;
    1.75 +    ValueVector _lower;
    1.76 +    ValueVector _upper;
    1.77 +    ValueVector _cap;
    1.78      CostVector _cost;
    1.79 -    FlowVector _supply;
    1.80 -    FlowVector _flow;
    1.81 +    ValueVector _supply;
    1.82 +    ValueVector _flow;
    1.83      CostVector _pi;
    1.84  
    1.85      // Data for storing the spanning tree structure
    1.86 @@ -689,12 +663,7 @@
    1.87      ///
    1.88      /// \param graph The digraph the algorithm runs on.
    1.89      NetworkSimplex(const GR& graph) :
    1.90 -      _graph(graph),
    1.91 -      _plower(NULL), _pupper(NULL), _pcost(NULL),
    1.92 -      _psupply(NULL), _pstsup(false), _stype(GEQ),
    1.93 -      _flow_map(NULL), _potential_map(NULL),
    1.94 -      _local_flow(false), _local_potential(false),
    1.95 -      _node_id(graph),
    1.96 +      _graph(graph), _node_id(graph), _arc_id(graph),
    1.97        INF(std::numeric_limits<Value>::has_infinity ?
    1.98            std::numeric_limits<Value>::infinity() :
    1.99            std::numeric_limits<Value>::max())
   1.100 @@ -704,12 +673,58 @@
   1.101          "The flow type of NetworkSimplex must be signed");
   1.102        LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   1.103          "The cost type of NetworkSimplex must be signed");
   1.104 -    }
   1.105 +        
   1.106 +      // Resize vectors
   1.107 +      _node_num = countNodes(_graph);
   1.108 +      _arc_num = countArcs(_graph);
   1.109 +      int all_node_num = _node_num + 1;
   1.110 +      int all_arc_num = _arc_num + _node_num;
   1.111  
   1.112 -    /// Destructor.
   1.113 -    ~NetworkSimplex() {
   1.114 -      if (_local_flow) delete _flow_map;
   1.115 -      if (_local_potential) delete _potential_map;
   1.116 +      _source.resize(all_arc_num);
   1.117 +      _target.resize(all_arc_num);
   1.118 +
   1.119 +      _lower.resize(all_arc_num);
   1.120 +      _upper.resize(all_arc_num);
   1.121 +      _cap.resize(all_arc_num);
   1.122 +      _cost.resize(all_arc_num);
   1.123 +      _supply.resize(all_node_num);
   1.124 +      _flow.resize(all_arc_num);
   1.125 +      _pi.resize(all_node_num);
   1.126 +
   1.127 +      _parent.resize(all_node_num);
   1.128 +      _pred.resize(all_node_num);
   1.129 +      _forward.resize(all_node_num);
   1.130 +      _thread.resize(all_node_num);
   1.131 +      _rev_thread.resize(all_node_num);
   1.132 +      _succ_num.resize(all_node_num);
   1.133 +      _last_succ.resize(all_node_num);
   1.134 +      _state.resize(all_arc_num);
   1.135 +
   1.136 +      // Copy the graph (store the arcs in a mixed order)
   1.137 +      int i = 0;
   1.138 +      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   1.139 +        _node_id[n] = i;
   1.140 +      }
   1.141 +      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
   1.142 +      i = 0;
   1.143 +      for (ArcIt a(_graph); a != INVALID; ++a) {
   1.144 +        _arc_id[a] = i;
   1.145 +        _source[i] = _node_id[_graph.source(a)];
   1.146 +        _target[i] = _node_id[_graph.target(a)];
   1.147 +        if ((i += k) >= _arc_num) i = (i % k) + 1;
   1.148 +      }
   1.149 +      
   1.150 +      // Initialize maps
   1.151 +      for (int i = 0; i != _node_num; ++i) {
   1.152 +        _supply[i] = 0;
   1.153 +      }
   1.154 +      for (int i = 0; i != _arc_num; ++i) {
   1.155 +        _lower[i] = 0;
   1.156 +        _upper[i] = INF;
   1.157 +        _cost[i] = 1;
   1.158 +      }
   1.159 +      _have_lower = false;
   1.160 +      _stype = GEQ;
   1.161      }
   1.162  
   1.163      /// \name Parameters
   1.164 @@ -731,10 +746,9 @@
   1.165      /// \return <tt>(*this)</tt>
   1.166      template <typename LowerMap>
   1.167      NetworkSimplex& lowerMap(const LowerMap& map) {
   1.168 -      delete _plower;
   1.169 -      _plower = new ValueArcMap(_graph);
   1.170 +      _have_lower = true;
   1.171        for (ArcIt a(_graph); a != INVALID; ++a) {
   1.172 -        (*_plower)[a] = map[a];
   1.173 +        _lower[_arc_id[a]] = map[a];
   1.174        }
   1.175        return *this;
   1.176      }
   1.177 @@ -753,10 +767,8 @@
   1.178      /// \return <tt>(*this)</tt>
   1.179      template<typename UpperMap>
   1.180      NetworkSimplex& upperMap(const UpperMap& map) {
   1.181 -      delete _pupper;
   1.182 -      _pupper = new ValueArcMap(_graph);
   1.183        for (ArcIt a(_graph); a != INVALID; ++a) {
   1.184 -        (*_pupper)[a] = map[a];
   1.185 +        _upper[_arc_id[a]] = map[a];
   1.186        }
   1.187        return *this;
   1.188      }
   1.189 @@ -774,10 +786,8 @@
   1.190      /// \return <tt>(*this)</tt>
   1.191      template<typename CostMap>
   1.192      NetworkSimplex& costMap(const CostMap& map) {
   1.193 -      delete _pcost;
   1.194 -      _pcost = new CostArcMap(_graph);
   1.195        for (ArcIt a(_graph); a != INVALID; ++a) {
   1.196 -        (*_pcost)[a] = map[a];
   1.197 +        _cost[_arc_id[a]] = map[a];
   1.198        }
   1.199        return *this;
   1.200      }
   1.201 @@ -796,11 +806,8 @@
   1.202      /// \return <tt>(*this)</tt>
   1.203      template<typename SupplyMap>
   1.204      NetworkSimplex& supplyMap(const SupplyMap& map) {
   1.205 -      delete _psupply;
   1.206 -      _pstsup = false;
   1.207 -      _psupply = new ValueNodeMap(_graph);
   1.208        for (NodeIt n(_graph); n != INVALID; ++n) {
   1.209 -        (*_psupply)[n] = map[n];
   1.210 +        _supply[_node_id[n]] = map[n];
   1.211        }
   1.212        return *this;
   1.213      }
   1.214 @@ -824,12 +831,11 @@
   1.215      ///
   1.216      /// \return <tt>(*this)</tt>
   1.217      NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
   1.218 -      delete _psupply;
   1.219 -      _psupply = NULL;
   1.220 -      _pstsup = true;
   1.221 -      _psource = s;
   1.222 -      _ptarget = t;
   1.223 -      _pstflow = k;
   1.224 +      for (int i = 0; i != _node_num; ++i) {
   1.225 +        _supply[i] = 0;
   1.226 +      }
   1.227 +      _supply[_node_id[s]] =  k;
   1.228 +      _supply[_node_id[t]] = -k;
   1.229        return *this;
   1.230      }
   1.231      
   1.232 @@ -847,41 +853,6 @@
   1.233        return *this;
   1.234      }
   1.235  
   1.236 -    /// \brief Set the flow map.
   1.237 -    ///
   1.238 -    /// This function sets the flow map.
   1.239 -    /// If it is not used before calling \ref run(), an instance will
   1.240 -    /// be allocated automatically. The destructor deallocates this
   1.241 -    /// automatically allocated map, of course.
   1.242 -    ///
   1.243 -    /// \return <tt>(*this)</tt>
   1.244 -    NetworkSimplex& flowMap(FlowMap& map) {
   1.245 -      if (_local_flow) {
   1.246 -        delete _flow_map;
   1.247 -        _local_flow = false;
   1.248 -      }
   1.249 -      _flow_map = &map;
   1.250 -      return *this;
   1.251 -    }
   1.252 -
   1.253 -    /// \brief Set the potential map.
   1.254 -    ///
   1.255 -    /// This function sets the potential map, which is used for storing
   1.256 -    /// the dual solution.
   1.257 -    /// If it is not used before calling \ref run(), an instance will
   1.258 -    /// be allocated automatically. The destructor deallocates this
   1.259 -    /// automatically allocated map, of course.
   1.260 -    ///
   1.261 -    /// \return <tt>(*this)</tt>
   1.262 -    NetworkSimplex& potentialMap(PotentialMap& map) {
   1.263 -      if (_local_potential) {
   1.264 -        delete _potential_map;
   1.265 -        _local_potential = false;
   1.266 -      }
   1.267 -      _potential_map = &map;
   1.268 -      return *this;
   1.269 -    }
   1.270 -    
   1.271      /// @}
   1.272  
   1.273      /// \name Execution Control
   1.274 @@ -894,7 +865,7 @@
   1.275      /// This function runs the algorithm.
   1.276      /// The paramters can be specified using functions \ref lowerMap(),
   1.277      /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   1.278 -    /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
   1.279 +    /// \ref supplyType().
   1.280      /// For example,
   1.281      /// \code
   1.282      ///   NetworkSimplex<ListDigraph> ns(graph);
   1.283 @@ -906,6 +877,8 @@
   1.284      /// that have been given are kept for the next call, unless
   1.285      /// \ref reset() is called, thus only the modified parameters
   1.286      /// have to be set again. See \ref reset() for examples.
   1.287 +    /// However the underlying digraph must not be modified after this
   1.288 +    /// class have been constructed, since it copies and extends the graph.
   1.289      ///
   1.290      /// \param pivot_rule The pivot rule that will be used during the
   1.291      /// algorithm. For more information see \ref PivotRule.
   1.292 @@ -928,12 +901,13 @@
   1.293      ///
   1.294      /// This function resets all the paramaters that have been given
   1.295      /// before using functions \ref lowerMap(), \ref upperMap(),
   1.296 -    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
   1.297 -    /// \ref flowMap() and \ref potentialMap().
   1.298 +    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
   1.299      ///
   1.300      /// It is useful for multiple run() calls. If this function is not
   1.301      /// used, all the parameters given before are kept for the next
   1.302      /// \ref run() call.
   1.303 +    /// However the underlying digraph must not be modified after this
   1.304 +    /// class have been constructed, since it copies and extends the graph.
   1.305      ///
   1.306      /// For example,
   1.307      /// \code
   1.308 @@ -957,23 +931,16 @@
   1.309      ///
   1.310      /// \return <tt>(*this)</tt>
   1.311      NetworkSimplex& reset() {
   1.312 -      delete _plower;
   1.313 -      delete _pupper;
   1.314 -      delete _pcost;
   1.315 -      delete _psupply;
   1.316 -      _plower = NULL;
   1.317 -      _pupper = NULL;
   1.318 -      _pcost = NULL;
   1.319 -      _psupply = NULL;
   1.320 -      _pstsup = false;
   1.321 +      for (int i = 0; i != _node_num; ++i) {
   1.322 +        _supply[i] = 0;
   1.323 +      }
   1.324 +      for (int i = 0; i != _arc_num; ++i) {
   1.325 +        _lower[i] = 0;
   1.326 +        _upper[i] = INF;
   1.327 +        _cost[i] = 1;
   1.328 +      }
   1.329 +      _have_lower = false;
   1.330        _stype = GEQ;
   1.331 -      if (_local_flow) delete _flow_map;
   1.332 -      if (_local_potential) delete _potential_map;
   1.333 -      _flow_map = NULL;
   1.334 -      _potential_map = NULL;
   1.335 -      _local_flow = false;
   1.336 -      _local_potential = false;
   1.337 -
   1.338        return *this;
   1.339      }
   1.340  
   1.341 @@ -1001,15 +968,12 @@
   1.342      /// function.
   1.343      ///
   1.344      /// \pre \ref run() must be called before using this function.
   1.345 -    template <typename Value>
   1.346 -    Value totalCost() const {
   1.347 -      Value c = 0;
   1.348 -      if (_pcost) {
   1.349 -        for (ArcIt e(_graph); e != INVALID; ++e)
   1.350 -          c += (*_flow_map)[e] * (*_pcost)[e];
   1.351 -      } else {
   1.352 -        for (ArcIt e(_graph); e != INVALID; ++e)
   1.353 -          c += (*_flow_map)[e];
   1.354 +    template <typename Number>
   1.355 +    Number totalCost() const {
   1.356 +      Number c = 0;
   1.357 +      for (ArcIt a(_graph); a != INVALID; ++a) {
   1.358 +        int i = _arc_id[a];
   1.359 +        c += Number(_flow[i]) * Number(_cost[i]);
   1.360        }
   1.361        return c;
   1.362      }
   1.363 @@ -1026,17 +990,21 @@
   1.364      ///
   1.365      /// \pre \ref run() must be called before using this function.
   1.366      Value flow(const Arc& a) const {
   1.367 -      return (*_flow_map)[a];
   1.368 +      return _flow[_arc_id[a]];
   1.369      }
   1.370  
   1.371 -    /// \brief Return a const reference to the flow map.
   1.372 +    /// \brief Return the flow map (the primal solution).
   1.373      ///
   1.374 -    /// This function returns a const reference to an arc map storing
   1.375 -    /// the found flow.
   1.376 +    /// This function copies the flow value on each arc into the given
   1.377 +    /// map. The \c Value type of the algorithm must be convertible to
   1.378 +    /// the \c Value type of the map.
   1.379      ///
   1.380      /// \pre \ref run() must be called before using this function.
   1.381 -    const FlowMap& flowMap() const {
   1.382 -      return *_flow_map;
   1.383 +    template <typename FlowMap>
   1.384 +    void flowMap(FlowMap &map) const {
   1.385 +      for (ArcIt a(_graph); a != INVALID; ++a) {
   1.386 +        map.set(a, _flow[_arc_id[a]]);
   1.387 +      }
   1.388      }
   1.389  
   1.390      /// \brief Return the potential (dual value) of the given node.
   1.391 @@ -1046,19 +1014,22 @@
   1.392      ///
   1.393      /// \pre \ref run() must be called before using this function.
   1.394      Cost potential(const Node& n) const {
   1.395 -      return (*_potential_map)[n];
   1.396 +      return _pi[_node_id[n]];
   1.397      }
   1.398  
   1.399 -    /// \brief Return a const reference to the potential map
   1.400 -    /// (the dual solution).
   1.401 +    /// \brief Return the potential map (the dual solution).
   1.402      ///
   1.403 -    /// This function returns a const reference to a node map storing
   1.404 -    /// the found potentials, which form the dual solution of the
   1.405 -    /// \ref min_cost_flow "minimum cost flow problem".
   1.406 +    /// This function copies the potential (dual value) of each node
   1.407 +    /// into the given map.
   1.408 +    /// The \c Cost type of the algorithm must be convertible to the
   1.409 +    /// \c Value type of the map.
   1.410      ///
   1.411      /// \pre \ref run() must be called before using this function.
   1.412 -    const PotentialMap& potentialMap() const {
   1.413 -      return *_potential_map;
   1.414 +    template <typename PotentialMap>
   1.415 +    void potentialMap(PotentialMap &map) const {
   1.416 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.417 +        map.set(n, _pi[_node_id[n]]);
   1.418 +      }
   1.419      }
   1.420  
   1.421      /// @}
   1.422 @@ -1067,69 +1038,33 @@
   1.423  
   1.424      // Initialize internal data structures
   1.425      bool init() {
   1.426 -      // Initialize result maps
   1.427 -      if (!_flow_map) {
   1.428 -        _flow_map = new FlowMap(_graph);
   1.429 -        _local_flow = true;
   1.430 -      }
   1.431 -      if (!_potential_map) {
   1.432 -        _potential_map = new PotentialMap(_graph);
   1.433 -        _local_potential = true;
   1.434 -      }
   1.435 -
   1.436 -      // Initialize vectors
   1.437 -      _node_num = countNodes(_graph);
   1.438 -      _arc_num = countArcs(_graph);
   1.439 -      int all_node_num = _node_num + 1;
   1.440 -      int all_arc_num = _arc_num + _node_num;
   1.441        if (_node_num == 0) return false;
   1.442  
   1.443 -      _arc_ref.resize(_arc_num);
   1.444 -      _source.resize(all_arc_num);
   1.445 -      _target.resize(all_arc_num);
   1.446 +      // Check the sum of supply values
   1.447 +      _sum_supply = 0;
   1.448 +      for (int i = 0; i != _node_num; ++i) {
   1.449 +        _sum_supply += _supply[i];
   1.450 +      }
   1.451 +      if ( !(_stype == GEQ && _sum_supply <= 0 ||
   1.452 +             _stype == LEQ && _sum_supply >= 0) ) return false;
   1.453  
   1.454 -      _cap.resize(all_arc_num);
   1.455 -      _cost.resize(all_arc_num);
   1.456 -      _supply.resize(all_node_num);
   1.457 -      _flow.resize(all_arc_num);
   1.458 -      _pi.resize(all_node_num);
   1.459 -
   1.460 -      _parent.resize(all_node_num);
   1.461 -      _pred.resize(all_node_num);
   1.462 -      _forward.resize(all_node_num);
   1.463 -      _thread.resize(all_node_num);
   1.464 -      _rev_thread.resize(all_node_num);
   1.465 -      _succ_num.resize(all_node_num);
   1.466 -      _last_succ.resize(all_node_num);
   1.467 -      _state.resize(all_arc_num);
   1.468 -
   1.469 -      // Initialize node related data
   1.470 -      bool valid_supply = true;
   1.471 -      _sum_supply = 0;
   1.472 -      if (!_pstsup && !_psupply) {
   1.473 -        _pstsup = true;
   1.474 -        _psource = _ptarget = NodeIt(_graph);
   1.475 -        _pstflow = 0;
   1.476 +      // Remove non-zero lower bounds
   1.477 +      if (_have_lower) {
   1.478 +        for (int i = 0; i != _arc_num; ++i) {
   1.479 +          Value c = _lower[i];
   1.480 +          if (c >= 0) {
   1.481 +            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
   1.482 +          } else {
   1.483 +            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
   1.484 +          }
   1.485 +          _supply[_source[i]] -= c;
   1.486 +          _supply[_target[i]] += c;
   1.487 +        }
   1.488 +      } else {
   1.489 +        for (int i = 0; i != _arc_num; ++i) {
   1.490 +          _cap[i] = _upper[i];
   1.491 +        }
   1.492        }
   1.493 -      if (_psupply) {
   1.494 -        int i = 0;
   1.495 -        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   1.496 -          _node_id[n] = i;
   1.497 -          _supply[i] = (*_psupply)[n];
   1.498 -          _sum_supply += _supply[i];
   1.499 -        }
   1.500 -        valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
   1.501 -                       (_stype == LEQ && _sum_supply >= 0);
   1.502 -      } else {
   1.503 -        int i = 0;
   1.504 -        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   1.505 -          _node_id[n] = i;
   1.506 -          _supply[i] = 0;
   1.507 -        }
   1.508 -        _supply[_node_id[_psource]] =  _pstflow;
   1.509 -        _supply[_node_id[_ptarget]] = -_pstflow;
   1.510 -      }
   1.511 -      if (!valid_supply) return false;
   1.512  
   1.513        // Initialize artifical cost
   1.514        Cost ART_COST;
   1.515 @@ -1143,93 +1078,31 @@
   1.516          ART_COST = (ART_COST + 1) * _node_num;
   1.517        }
   1.518  
   1.519 +      // Initialize arc maps
   1.520 +      for (int i = 0; i != _arc_num; ++i) {
   1.521 +        _flow[i] = 0;
   1.522 +        _state[i] = STATE_LOWER;
   1.523 +      }
   1.524 +      
   1.525        // Set data for the artificial root node
   1.526        _root = _node_num;
   1.527        _parent[_root] = -1;
   1.528        _pred[_root] = -1;
   1.529        _thread[_root] = 0;
   1.530        _rev_thread[0] = _root;
   1.531 -      _succ_num[_root] = all_node_num;
   1.532 +      _succ_num[_root] = _node_num + 1;
   1.533        _last_succ[_root] = _root - 1;
   1.534        _supply[_root] = -_sum_supply;
   1.535 -      if (_sum_supply < 0) {
   1.536 -        _pi[_root] = -ART_COST;
   1.537 -      } else {
   1.538 -        _pi[_root] = ART_COST;
   1.539 -      }
   1.540 -
   1.541 -      // Store the arcs in a mixed order
   1.542 -      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
   1.543 -      int i = 0;
   1.544 -      for (ArcIt e(_graph); e != INVALID; ++e) {
   1.545 -        _arc_ref[i] = e;
   1.546 -        if ((i += k) >= _arc_num) i = (i % k) + 1;
   1.547 -      }
   1.548 -
   1.549 -      // Initialize arc maps
   1.550 -      if (_pupper && _pcost) {
   1.551 -        for (int i = 0; i != _arc_num; ++i) {
   1.552 -          Arc e = _arc_ref[i];
   1.553 -          _source[i] = _node_id[_graph.source(e)];
   1.554 -          _target[i] = _node_id[_graph.target(e)];
   1.555 -          _cap[i] = (*_pupper)[e];
   1.556 -          _cost[i] = (*_pcost)[e];
   1.557 -          _flow[i] = 0;
   1.558 -          _state[i] = STATE_LOWER;
   1.559 -        }
   1.560 -      } else {
   1.561 -        for (int i = 0; i != _arc_num; ++i) {
   1.562 -          Arc e = _arc_ref[i];
   1.563 -          _source[i] = _node_id[_graph.source(e)];
   1.564 -          _target[i] = _node_id[_graph.target(e)];
   1.565 -          _flow[i] = 0;
   1.566 -          _state[i] = STATE_LOWER;
   1.567 -        }
   1.568 -        if (_pupper) {
   1.569 -          for (int i = 0; i != _arc_num; ++i)
   1.570 -            _cap[i] = (*_pupper)[_arc_ref[i]];
   1.571 -        } else {
   1.572 -          for (int i = 0; i != _arc_num; ++i)
   1.573 -            _cap[i] = INF;
   1.574 -        }
   1.575 -        if (_pcost) {
   1.576 -          for (int i = 0; i != _arc_num; ++i)
   1.577 -            _cost[i] = (*_pcost)[_arc_ref[i]];
   1.578 -        } else {
   1.579 -          for (int i = 0; i != _arc_num; ++i)
   1.580 -            _cost[i] = 1;
   1.581 -        }
   1.582 -      }
   1.583 -      
   1.584 -      // Remove non-zero lower bounds
   1.585 -      if (_plower) {
   1.586 -        for (int i = 0; i != _arc_num; ++i) {
   1.587 -          Value c = (*_plower)[_arc_ref[i]];
   1.588 -          if (c > 0) {
   1.589 -            if (_cap[i] < INF) _cap[i] -= c;
   1.590 -            _supply[_source[i]] -= c;
   1.591 -            _supply[_target[i]] += c;
   1.592 -          }
   1.593 -          else if (c < 0) {
   1.594 -            if (_cap[i] < INF + c) {
   1.595 -              _cap[i] -= c;
   1.596 -            } else {
   1.597 -              _cap[i] = INF;
   1.598 -            }
   1.599 -            _supply[_source[i]] -= c;
   1.600 -            _supply[_target[i]] += c;
   1.601 -          }
   1.602 -        }
   1.603 -      }
   1.604 +      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
   1.605  
   1.606        // Add artificial arcs and initialize the spanning tree data structure
   1.607        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   1.608 +        _parent[u] = _root;
   1.609 +        _pred[u] = e;
   1.610          _thread[u] = u + 1;
   1.611          _rev_thread[u + 1] = u;
   1.612          _succ_num[u] = 1;
   1.613          _last_succ[u] = u;
   1.614 -        _parent[u] = _root;
   1.615 -        _pred[u] = e;
   1.616          _cost[e] = ART_COST;
   1.617          _cap[e] = INF;
   1.618          _state[e] = STATE_TREE;
   1.619 @@ -1517,20 +1390,16 @@
   1.620          }
   1.621        }
   1.622  
   1.623 -      // Copy flow values to _flow_map
   1.624 -      if (_plower) {
   1.625 +      // Transform the solution and the supply map to the original form
   1.626 +      if (_have_lower) {
   1.627          for (int i = 0; i != _arc_num; ++i) {
   1.628 -          Arc e = _arc_ref[i];
   1.629 -          _flow_map->set(e, (*_plower)[e] + _flow[i]);
   1.630 +          Value c = _lower[i];
   1.631 +          if (c != 0) {
   1.632 +            _flow[i] += c;
   1.633 +            _supply[_source[i]] += c;
   1.634 +            _supply[_target[i]] -= c;
   1.635 +          }
   1.636          }
   1.637 -      } else {
   1.638 -        for (int i = 0; i != _arc_num; ++i) {
   1.639 -          _flow_map->set(_arc_ref[i], _flow[i]);
   1.640 -        }
   1.641 -      }
   1.642 -      // Copy potential values to _potential_map
   1.643 -      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.644 -        _potential_map->set(n, _pi[_node_id[n]]);
   1.645        }
   1.646  
   1.647        return OPTIMAL;
     2.1 --- a/test/min_cost_flow_test.cc	Wed Apr 29 14:25:51 2009 +0200
     2.2 +++ b/test/min_cost_flow_test.cc	Wed Apr 29 16:54:27 2009 +0200
     2.3 @@ -84,7 +84,7 @@
     2.4  };
     2.5  
     2.6  // Check the interface of an MCF algorithm
     2.7 -template <typename GR, typename Flow, typename Cost>
     2.8 +template <typename GR, typename Value, typename Cost>
     2.9  class McfClassConcept
    2.10  {
    2.11  public:
    2.12 @@ -95,6 +95,7 @@
    2.13        checkConcept<concepts::Digraph, GR>();
    2.14  
    2.15        MCF mcf(g);
    2.16 +      const MCF& const_mcf = mcf;
    2.17  
    2.18        b = mcf.reset()
    2.19               .lowerMap(lower)
    2.20 @@ -102,47 +103,38 @@
    2.21               .costMap(cost)
    2.22               .supplyMap(sup)
    2.23               .stSupply(n, n, k)
    2.24 -             .flowMap(flow)
    2.25 -             .potentialMap(pot)
    2.26               .run();
    2.27 -      
    2.28 -      const MCF& const_mcf = mcf;
    2.29 -
    2.30 -      const typename MCF::FlowMap &fm = const_mcf.flowMap();
    2.31 -      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
    2.32  
    2.33        c = const_mcf.totalCost();
    2.34 -      double x = const_mcf.template totalCost<double>();
    2.35 +      x = const_mcf.template totalCost<double>();
    2.36        v = const_mcf.flow(a);
    2.37        c = const_mcf.potential(n);
    2.38 -      
    2.39 -      v = const_mcf.INF;
    2.40 -
    2.41 -      ignore_unused_variable_warning(fm);
    2.42 -      ignore_unused_variable_warning(pm);
    2.43 -      ignore_unused_variable_warning(x);
    2.44 +      const_mcf.flowMap(fm);
    2.45 +      const_mcf.potentialMap(pm);
    2.46      }
    2.47  
    2.48      typedef typename GR::Node Node;
    2.49      typedef typename GR::Arc Arc;
    2.50 -    typedef concepts::ReadMap<Node, Flow> NM;
    2.51 -    typedef concepts::ReadMap<Arc, Flow> FAM;
    2.52 +    typedef concepts::ReadMap<Node, Value> NM;
    2.53 +    typedef concepts::ReadMap<Arc, Value> VAM;
    2.54      typedef concepts::ReadMap<Arc, Cost> CAM;
    2.55 +    typedef concepts::WriteMap<Arc, Value> FlowMap;
    2.56 +    typedef concepts::WriteMap<Node, Cost> PotMap;
    2.57  
    2.58      const GR &g;
    2.59 -    const FAM &lower;
    2.60 -    const FAM &upper;
    2.61 +    const VAM &lower;
    2.62 +    const VAM &upper;
    2.63      const CAM &cost;
    2.64      const NM &sup;
    2.65      const Node &n;
    2.66      const Arc &a;
    2.67 -    const Flow &k;
    2.68 -    Flow v;
    2.69 -    Cost c;
    2.70 +    const Value &k;
    2.71 +    FlowMap fm;
    2.72 +    PotMap pm;
    2.73      bool b;
    2.74 -
    2.75 -    typename MCF::FlowMap &flow;
    2.76 -    typename MCF::PotentialMap &pot;
    2.77 +    double x;
    2.78 +    typename MCF::Value v;
    2.79 +    typename MCF::Cost c;
    2.80    };
    2.81  
    2.82  };
    2.83 @@ -221,11 +213,14 @@
    2.84  {
    2.85    check(mcf_result == result, "Wrong result " + test_id);
    2.86    if (optimal) {
    2.87 -    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
    2.88 +    typename GR::template ArcMap<typename SM::Value> flow(gr);
    2.89 +    typename GR::template NodeMap<typename CM::Value> pi(gr);
    2.90 +    mcf.flowMap(flow);
    2.91 +    mcf.potentialMap(pi);
    2.92 +    check(checkFlow(gr, lower, upper, supply, flow, type),
    2.93            "The flow is not feasible " + test_id);
    2.94      check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
    2.95 -    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
    2.96 -                         mcf.potentialMap()),
    2.97 +    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
    2.98            "Wrong potentials " + test_id);
    2.99    }
   2.100  }
   2.101 @@ -234,11 +229,13 @@
   2.102  {
   2.103    // Check the interfaces
   2.104    {
   2.105 -    typedef int Flow;
   2.106 -    typedef int Cost;
   2.107      typedef concepts::Digraph GR;
   2.108 -    checkConcept< McfClassConcept<GR, Flow, Cost>,
   2.109 -                  NetworkSimplex<GR, Flow, Cost> >();
   2.110 +    checkConcept< McfClassConcept<GR, int, int>,
   2.111 +                  NetworkSimplex<GR> >();
   2.112 +    checkConcept< McfClassConcept<GR, double, double>,
   2.113 +                  NetworkSimplex<GR, double> >();
   2.114 +    checkConcept< McfClassConcept<GR, int, double>,
   2.115 +                  NetworkSimplex<GR, int, double> >();
   2.116    }
   2.117  
   2.118    // Run various MCF tests