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 = ↦
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 = ↦
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 ⊃
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