Query improvements in the min cost flow algorithms.
- External flow and potential maps can be used.
- New query functions: flow() and potential().
- CycleCanceling also provides dual solution (node potentials).
- Doc improvements.
1.1 --- a/lemon/capacity_scaling.h Wed Feb 27 11:39:03 2008 +0000
1.2 +++ b/lemon/capacity_scaling.h Thu Feb 28 02:54:27 2008 +0000
1.3 @@ -50,11 +50,8 @@
1.4 /// \warning
1.5 /// - Edge capacities and costs should be \e non-negative \e integers.
1.6 /// - Supply values should be \e signed \e integers.
1.7 - /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
1.8 - /// - \c CapacityMap::Value and \c SupplyMap::Value must be
1.9 - /// convertible to each other.
1.10 - /// - All value types must be convertible to \c CostMap::Value, which
1.11 - /// must be signed type.
1.12 + /// - The value types of the maps should be convertible to each other.
1.13 + /// - \c CostMap::Value must be signed type.
1.14 ///
1.15 /// \author Peter Kovacs
1.16
1.17 @@ -120,7 +117,7 @@
1.18
1.19 public:
1.20
1.21 - /// The constructor of the class.
1.22 + /// Constructor.
1.23 ResidualDijkstra( const Graph &graph,
1.24 const FlowMap &flow,
1.25 const CapacityEdgeMap &res_cap,
1.26 @@ -221,9 +218,11 @@
1.27 bool _valid_supply;
1.28
1.29 // Edge map of the current flow
1.30 - FlowMap _flow;
1.31 + FlowMap *_flow;
1.32 + bool _local_flow;
1.33 // Node map of the current potentials
1.34 - PotentialMap _potential;
1.35 + PotentialMap *_potential;
1.36 + bool _local_potential;
1.37
1.38 // The residual capacity map
1.39 CapacityEdgeMap _res_cap;
1.40 @@ -243,13 +242,13 @@
1.41 PredMap _pred;
1.42 // Implementation of the Dijkstra algorithm for finding augmenting
1.43 // shortest paths in the residual network
1.44 - ResidualDijkstra _dijkstra;
1.45 + ResidualDijkstra *_dijkstra;
1.46
1.47 - public :
1.48 + public:
1.49
1.50 - /// \brief General constructor of the class (with lower bounds).
1.51 + /// \brief General constructor (with lower bounds).
1.52 ///
1.53 - /// General constructor of the class (with lower bounds).
1.54 + /// General constructor (with lower bounds).
1.55 ///
1.56 /// \param graph The directed graph the algorithm runs on.
1.57 /// \param lower The lower bounds of the edges.
1.58 @@ -262,9 +261,9 @@
1.59 const CostMap &cost,
1.60 const SupplyMap &supply ) :
1.61 _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
1.62 - _supply(graph), _flow(graph, 0), _potential(graph, 0),
1.63 - _res_cap(graph), _excess(graph), _pred(graph),
1.64 - _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
1.65 + _supply(graph), _flow(0), _local_flow(false),
1.66 + _potential(0), _local_potential(false),
1.67 + _res_cap(graph), _excess(graph), _pred(graph)
1.68 {
1.69 // Removing non-zero lower bounds
1.70 _capacity = subMap(capacity, lower);
1.71 @@ -282,9 +281,9 @@
1.72 _valid_supply = sum == 0;
1.73 }
1.74
1.75 - /// \brief General constructor of the class (without lower bounds).
1.76 + /// \brief General constructor (without lower bounds).
1.77 ///
1.78 - /// General constructor of the class (without lower bounds).
1.79 + /// General constructor (without lower bounds).
1.80 ///
1.81 /// \param graph The directed graph the algorithm runs on.
1.82 /// \param capacity The capacities (upper bounds) of the edges.
1.83 @@ -295,9 +294,9 @@
1.84 const CostMap &cost,
1.85 const SupplyMap &supply ) :
1.86 _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
1.87 - _supply(supply), _flow(graph, 0), _potential(graph, 0),
1.88 - _res_cap(capacity), _excess(graph), _pred(graph),
1.89 - _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
1.90 + _supply(supply), _flow(0), _local_flow(false),
1.91 + _potential(0), _local_potential(false),
1.92 + _res_cap(capacity), _excess(graph), _pred(graph)
1.93 {
1.94 // Checking the sum of supply values
1.95 Supply sum = 0;
1.96 @@ -305,9 +304,9 @@
1.97 _valid_supply = sum == 0;
1.98 }
1.99
1.100 - /// \brief Simple constructor of the class (with lower bounds).
1.101 + /// \brief Simple constructor (with lower bounds).
1.102 ///
1.103 - /// Simple constructor of the class (with lower bounds).
1.104 + /// Simple constructor (with lower bounds).
1.105 ///
1.106 /// \param graph The directed graph the algorithm runs on.
1.107 /// \param lower The lower bounds of the edges.
1.108 @@ -324,9 +323,9 @@
1.109 Node s, Node t,
1.110 Supply flow_value ) :
1.111 _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
1.112 - _supply(graph), _flow(graph, 0), _potential(graph, 0),
1.113 - _res_cap(graph), _excess(graph), _pred(graph),
1.114 - _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
1.115 + _supply(graph), _flow(0), _local_flow(false),
1.116 + _potential(0), _local_potential(false),
1.117 + _res_cap(graph), _excess(graph), _pred(graph)
1.118 {
1.119 // Removing non-zero lower bounds
1.120 _capacity = subMap(capacity, lower);
1.121 @@ -344,9 +343,9 @@
1.122 _valid_supply = true;
1.123 }
1.124
1.125 - /// \brief Simple constructor of the class (without lower bounds).
1.126 + /// \brief Simple constructor (without lower bounds).
1.127 ///
1.128 - /// Simple constructor of the class (without lower bounds).
1.129 + /// Simple constructor (without lower bounds).
1.130 ///
1.131 /// \param graph The directed graph the algorithm runs on.
1.132 /// \param capacity The capacities (upper bounds) of the edges.
1.133 @@ -361,15 +360,56 @@
1.134 Node s, Node t,
1.135 Supply flow_value ) :
1.136 _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
1.137 - _supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
1.138 - _res_cap(capacity), _excess(graph), _pred(graph),
1.139 - _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
1.140 + _supply(graph, 0), _flow(0), _local_flow(false),
1.141 + _potential(0), _local_potential(false),
1.142 + _res_cap(capacity), _excess(graph), _pred(graph)
1.143 {
1.144 _supply[s] = flow_value;
1.145 _supply[t] = -flow_value;
1.146 _valid_supply = true;
1.147 }
1.148
1.149 + /// Destructor.
1.150 + ~CapacityScaling() {
1.151 + if (_local_flow) delete _flow;
1.152 + if (_local_potential) delete _potential;
1.153 + delete _dijkstra;
1.154 + }
1.155 +
1.156 + /// \brief Sets the flow map.
1.157 + ///
1.158 + /// Sets the flow map.
1.159 + ///
1.160 + /// \return \c (*this)
1.161 + CapacityScaling& flowMap(FlowMap &map) {
1.162 + if (_local_flow) {
1.163 + delete _flow;
1.164 + _local_flow = false;
1.165 + }
1.166 + _flow = ↦
1.167 + return *this;
1.168 + }
1.169 +
1.170 + /// \brief Sets the potential map.
1.171 + ///
1.172 + /// Sets the potential map.
1.173 + ///
1.174 + /// \return \c (*this)
1.175 + CapacityScaling& potentialMap(PotentialMap &map) {
1.176 + if (_local_potential) {
1.177 + delete _potential;
1.178 + _local_potential = false;
1.179 + }
1.180 + _potential = ↦
1.181 + return *this;
1.182 + }
1.183 +
1.184 + /// \name Execution control
1.185 + /// The only way to execute the algorithm is to call the run()
1.186 + /// function.
1.187 +
1.188 + /// @{
1.189 +
1.190 /// \brief Runs the algorithm.
1.191 ///
1.192 /// Runs the algorithm.
1.193 @@ -384,6 +424,15 @@
1.194 return init(scaling) && start();
1.195 }
1.196
1.197 + /// @}
1.198 +
1.199 + /// \name Query Functions
1.200 + /// The result of the algorithm can be obtained using these
1.201 + /// functions.
1.202 + /// \n run() must be called before using them.
1.203 +
1.204 + /// @{
1.205 +
1.206 /// \brief Returns a const reference to the edge map storing the
1.207 /// found flow.
1.208 ///
1.209 @@ -391,7 +440,7 @@
1.210 ///
1.211 /// \pre \ref run() must be called before using this function.
1.212 const FlowMap& flowMap() const {
1.213 - return _flow;
1.214 + return *_flow;
1.215 }
1.216
1.217 /// \brief Returns a const reference to the node map storing the
1.218 @@ -402,7 +451,25 @@
1.219 ///
1.220 /// \pre \ref run() must be called before using this function.
1.221 const PotentialMap& potentialMap() const {
1.222 - return _potential;
1.223 + return *_potential;
1.224 + }
1.225 +
1.226 + /// \brief Returns the flow on the edge.
1.227 + ///
1.228 + /// Returns the flow on the edge.
1.229 + ///
1.230 + /// \pre \ref run() must be called before using this function.
1.231 + Capacity flow(const Edge& edge) const {
1.232 + return (*_flow)[edge];
1.233 + }
1.234 +
1.235 + /// \brief Returns the potential of the node.
1.236 + ///
1.237 + /// Returns the potential of the node.
1.238 + ///
1.239 + /// \pre \ref run() must be called before using this function.
1.240 + Cost potential(const Node& node) const {
1.241 + return (*_potential)[node];
1.242 }
1.243
1.244 /// \brief Returns the total cost of the found flow.
1.245 @@ -414,18 +481,35 @@
1.246 Cost totalCost() const {
1.247 Cost c = 0;
1.248 for (EdgeIt e(_graph); e != INVALID; ++e)
1.249 - c += _flow[e] * _cost[e];
1.250 + c += (*_flow)[e] * _cost[e];
1.251 return c;
1.252 }
1.253
1.254 + /// @}
1.255 +
1.256 private:
1.257
1.258 /// Initializes the algorithm.
1.259 bool init(bool scaling) {
1.260 if (!_valid_supply) return false;
1.261 +
1.262 + // Initializing maps
1.263 + if (!_flow) {
1.264 + _flow = new FlowMap(_graph);
1.265 + _local_flow = true;
1.266 + }
1.267 + if (!_potential) {
1.268 + _potential = new PotentialMap(_graph);
1.269 + _local_potential = true;
1.270 + }
1.271 + for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
1.272 + for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
1.273 _excess = _supply;
1.274
1.275 - // Initilaizing delta value
1.276 + _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
1.277 + _excess, *_potential, _pred );
1.278 +
1.279 + // Initializing delta value
1.280 if (scaling) {
1.281 // With scaling
1.282 Supply max_sup = 0, max_dem = 0;
1.283 @@ -441,6 +525,7 @@
1.284 // Without scaling
1.285 _delta = 1;
1.286 }
1.287 +
1.288 return true;
1.289 }
1.290
1.291 @@ -461,17 +546,17 @@
1.292 // Saturating all edges not satisfying the optimality condition
1.293 for (EdgeIt e(_graph); e != INVALID; ++e) {
1.294 Node u = _graph.source(e), v = _graph.target(e);
1.295 - Cost c = _cost[e] + _potential[u] - _potential[v];
1.296 + Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
1.297 if (c < 0 && _res_cap[e] >= _delta) {
1.298 _excess[u] -= _res_cap[e];
1.299 _excess[v] += _res_cap[e];
1.300 - _flow[e] = _capacity[e];
1.301 + (*_flow)[e] = _capacity[e];
1.302 _res_cap[e] = 0;
1.303 }
1.304 - else if (c > 0 && _flow[e] >= _delta) {
1.305 - _excess[u] += _flow[e];
1.306 - _excess[v] -= _flow[e];
1.307 - _flow[e] = 0;
1.308 + else if (c > 0 && (*_flow)[e] >= _delta) {
1.309 + _excess[u] += (*_flow)[e];
1.310 + _excess[v] -= (*_flow)[e];
1.311 + (*_flow)[e] = 0;
1.312 _res_cap[e] = _capacity[e];
1.313 }
1.314 }
1.315 @@ -501,7 +586,7 @@
1.316
1.317 // Running Dijkstra
1.318 s = _excess_nodes[next_node];
1.319 - if ((t = _dijkstra.run(s, _delta)) == INVALID) {
1.320 + if ((t = _dijkstra->run(s, _delta)) == INVALID) {
1.321 if (_delta > 1) {
1.322 ++next_node;
1.323 continue;
1.324 @@ -520,7 +605,7 @@
1.325 rc = _res_cap[e];
1.326 u = _graph.source(e);
1.327 } else {
1.328 - rc = _flow[e];
1.329 + rc = (*_flow)[e];
1.330 u = _graph.target(e);
1.331 }
1.332 if (rc < d) d = rc;
1.333 @@ -529,11 +614,11 @@
1.334 u = t;
1.335 while ((e = _pred[u]) != INVALID) {
1.336 if (u == _graph.target(e)) {
1.337 - _flow[e] += d;
1.338 + (*_flow)[e] += d;
1.339 _res_cap[e] -= d;
1.340 u = _graph.source(e);
1.341 } else {
1.342 - _flow[e] -= d;
1.343 + (*_flow)[e] -= d;
1.344 _res_cap[e] += d;
1.345 u = _graph.target(e);
1.346 }
1.347 @@ -552,7 +637,7 @@
1.348 // Handling non-zero lower bounds
1.349 if (_lower) {
1.350 for (EdgeIt e(_graph); e != INVALID; ++e)
1.351 - _flow[e] += (*_lower)[e];
1.352 + (*_flow)[e] += (*_lower)[e];
1.353 }
1.354 return true;
1.355 }
1.356 @@ -572,7 +657,7 @@
1.357 {
1.358 // Running Dijkstra
1.359 s = _excess_nodes[next_node];
1.360 - if ((t = _dijkstra.run(s, 1)) == INVALID)
1.361 + if ((t = _dijkstra->run(s, 1)) == INVALID)
1.362 return false;
1.363
1.364 // Augmenting along a shortest path from s to t
1.365 @@ -585,7 +670,7 @@
1.366 rc = _res_cap[e];
1.367 u = _graph.source(e);
1.368 } else {
1.369 - rc = _flow[e];
1.370 + rc = (*_flow)[e];
1.371 u = _graph.target(e);
1.372 }
1.373 if (rc < d) d = rc;
1.374 @@ -593,11 +678,11 @@
1.375 u = t;
1.376 while ((e = _pred[u]) != INVALID) {
1.377 if (u == _graph.target(e)) {
1.378 - _flow[e] += d;
1.379 + (*_flow)[e] += d;
1.380 _res_cap[e] -= d;
1.381 u = _graph.source(e);
1.382 } else {
1.383 - _flow[e] -= d;
1.384 + (*_flow)[e] -= d;
1.385 _res_cap[e] += d;
1.386 u = _graph.target(e);
1.387 }
1.388 @@ -609,7 +694,7 @@
1.389 // Handling non-zero lower bounds
1.390 if (_lower) {
1.391 for (EdgeIt e(_graph); e != INVALID; ++e)
1.392 - _flow[e] += (*_lower)[e];
1.393 + (*_flow)[e] += (*_lower)[e];
1.394 }
1.395 return true;
1.396 }
2.1 --- a/lemon/cost_scaling.h Wed Feb 27 11:39:03 2008 +0000
2.2 +++ b/lemon/cost_scaling.h Thu Feb 28 02:54:27 2008 +0000
2.3 @@ -54,11 +54,8 @@
2.4 /// \warning
2.5 /// - Edge capacities and costs should be \e non-negative \e integers.
2.6 /// - Supply values should be \e signed \e integers.
2.7 - /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
2.8 - /// - \c CapacityMap::Value and \c SupplyMap::Value must be
2.9 - /// convertible to each other.
2.10 - /// - All value types must be convertible to \c CostMap::Value, which
2.11 - /// must be signed type.
2.12 + /// - The value types of the maps should be convertible to each other.
2.13 + /// - \c CostMap::Value must be signed type.
2.14 ///
2.15 /// \note Edge costs are multiplied with the number of nodes during
2.16 /// the algorithm so overflow problems may arise more easily than with
2.17 @@ -97,7 +94,7 @@
2.18 public:
2.19
2.20 /// The type of the flow map.
2.21 - typedef CapacityEdgeMap FlowMap;
2.22 + typedef typename Graph::template EdgeMap<Capacity> FlowMap;
2.23 /// The type of the potential map.
2.24 typedef typename Graph::template NodeMap<LCost> PotentialMap;
2.25
2.26 @@ -107,20 +104,21 @@
2.27 ///
2.28 /// \ref ResidualCostMap is a map adaptor class for handling
2.29 /// residual edge costs.
2.30 - class ResidualCostMap : public MapBase<ResEdge, LCost>
2.31 + template <typename Map>
2.32 + class ResidualCostMap : public MapBase<ResEdge, typename Map::Value>
2.33 {
2.34 private:
2.35
2.36 - const LargeCostMap &_cost_map;
2.37 + const Map &_cost_map;
2.38
2.39 public:
2.40
2.41 ///\e
2.42 - ResidualCostMap(const LargeCostMap &cost_map) :
2.43 + ResidualCostMap(const Map &cost_map) :
2.44 _cost_map(cost_map) {}
2.45
2.46 ///\e
2.47 - LCost operator[](const ResEdge &e) const {
2.48 + typename Map::Value operator[](const ResEdge &e) const {
2.49 return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
2.50 }
2.51
2.52 @@ -160,8 +158,8 @@
2.53 static const int ALPHA = 4;
2.54
2.55 // Paramters for heuristics
2.56 - static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
2.57 - static const int BF_HEURISTIC_BOUND_FACTOR = 3;
2.58 + static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
2.59 + static const int BF_HEURISTIC_BOUND_FACTOR = 3;
2.60
2.61 private:
2.62
2.63 @@ -180,16 +178,18 @@
2.64 bool _valid_supply;
2.65
2.66 // Edge map of the current flow
2.67 - FlowMap _flow;
2.68 + FlowMap *_flow;
2.69 + bool _local_flow;
2.70 // Node map of the current potentials
2.71 - PotentialMap _potential;
2.72 + PotentialMap *_potential;
2.73 + bool _local_potential;
2.74
2.75 // The residual graph
2.76 - ResGraph _res_graph;
2.77 + ResGraph *_res_graph;
2.78 // The residual cost map
2.79 - ResidualCostMap _res_cost;
2.80 + ResidualCostMap<LargeCostMap> _res_cost;
2.81 // The reduced cost map
2.82 - ReducedCostMap _red_cost;
2.83 + ReducedCostMap *_red_cost;
2.84 // The excess map
2.85 SupplyNodeMap _excess;
2.86 // The epsilon parameter used for cost scaling
2.87 @@ -197,9 +197,9 @@
2.88
2.89 public:
2.90
2.91 - /// \brief General constructor of the class (with lower bounds).
2.92 + /// \brief General constructor (with lower bounds).
2.93 ///
2.94 - /// General constructor of the class (with lower bounds).
2.95 + /// General constructor (with lower bounds).
2.96 ///
2.97 /// \param graph The directed graph the algorithm runs on.
2.98 /// \param lower The lower bounds of the edges.
2.99 @@ -212,9 +212,9 @@
2.100 const CostMap &cost,
2.101 const SupplyMap &supply ) :
2.102 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
2.103 - _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
2.104 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
2.105 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
2.106 + _cost(graph), _supply(graph), _flow(0), _local_flow(false),
2.107 + _potential(0), _local_potential(false), _res_cost(_cost),
2.108 + _excess(graph, 0)
2.109 {
2.110 // Removing non-zero lower bounds
2.111 _capacity = subMap(capacity, lower);
2.112 @@ -231,9 +231,9 @@
2.113 _valid_supply = sum == 0;
2.114 }
2.115
2.116 - /// \brief General constructor of the class (without lower bounds).
2.117 + /// \brief General constructor (without lower bounds).
2.118 ///
2.119 - /// General constructor of the class (without lower bounds).
2.120 + /// General constructor (without lower bounds).
2.121 ///
2.122 /// \param graph The directed graph the algorithm runs on.
2.123 /// \param capacity The capacities (upper bounds) of the edges.
2.124 @@ -244,9 +244,9 @@
2.125 const CostMap &cost,
2.126 const SupplyMap &supply ) :
2.127 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
2.128 - _cost(graph), _supply(supply), _flow(graph, 0), _potential(graph, 0),
2.129 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
2.130 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
2.131 + _cost(graph), _supply(supply), _flow(0), _local_flow(false),
2.132 + _potential(0), _local_potential(false), _res_cost(_cost),
2.133 + _excess(graph, 0)
2.134 {
2.135 // Checking the sum of supply values
2.136 Supply sum = 0;
2.137 @@ -254,9 +254,9 @@
2.138 _valid_supply = sum == 0;
2.139 }
2.140
2.141 - /// \brief Simple constructor of the class (with lower bounds).
2.142 + /// \brief Simple constructor (with lower bounds).
2.143 ///
2.144 - /// Simple constructor of the class (with lower bounds).
2.145 + /// Simple constructor (with lower bounds).
2.146 ///
2.147 /// \param graph The directed graph the algorithm runs on.
2.148 /// \param lower The lower bounds of the edges.
2.149 @@ -273,9 +273,9 @@
2.150 Node s, Node t,
2.151 Supply flow_value ) :
2.152 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
2.153 - _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
2.154 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
2.155 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
2.156 + _cost(graph), _supply(graph), _flow(0), _local_flow(false),
2.157 + _potential(0), _local_potential(false), _res_cost(_cost),
2.158 + _excess(graph, 0)
2.159 {
2.160 // Removing nonzero lower bounds
2.161 _capacity = subMap(capacity, lower);
2.162 @@ -292,9 +292,9 @@
2.163 _valid_supply = true;
2.164 }
2.165
2.166 - /// \brief Simple constructor of the class (without lower bounds).
2.167 + /// \brief Simple constructor (without lower bounds).
2.168 ///
2.169 - /// Simple constructor of the class (without lower bounds).
2.170 + /// Simple constructor (without lower bounds).
2.171 ///
2.172 /// \param graph The directed graph the algorithm runs on.
2.173 /// \param capacity The capacities (upper bounds) of the edges.
2.174 @@ -309,24 +309,75 @@
2.175 Node s, Node t,
2.176 Supply flow_value ) :
2.177 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
2.178 - _cost(graph), _supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
2.179 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
2.180 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
2.181 + _cost(graph), _supply(graph, 0), _flow(0), _local_flow(false),
2.182 + _potential(0), _local_potential(false), _res_cost(_cost),
2.183 + _excess(graph, 0)
2.184 {
2.185 _supply[s] = flow_value;
2.186 _supply[t] = -flow_value;
2.187 _valid_supply = true;
2.188 }
2.189
2.190 + /// Destructor.
2.191 + ~CostScaling() {
2.192 + if (_local_flow) delete _flow;
2.193 + if (_local_potential) delete _potential;
2.194 + delete _res_graph;
2.195 + delete _red_cost;
2.196 + }
2.197 +
2.198 + /// \brief Sets the flow map.
2.199 + ///
2.200 + /// Sets the flow map.
2.201 + ///
2.202 + /// \return \c (*this)
2.203 + CostScaling& flowMap(FlowMap &map) {
2.204 + if (_local_flow) {
2.205 + delete _flow;
2.206 + _local_flow = false;
2.207 + }
2.208 + _flow = ↦
2.209 + return *this;
2.210 + }
2.211 +
2.212 + /// \brief Sets the potential map.
2.213 + ///
2.214 + /// Sets the potential map.
2.215 + ///
2.216 + /// \return \c (*this)
2.217 + CostScaling& potentialMap(PotentialMap &map) {
2.218 + if (_local_potential) {
2.219 + delete _potential;
2.220 + _local_potential = false;
2.221 + }
2.222 + _potential = ↦
2.223 + return *this;
2.224 + }
2.225 +
2.226 + /// \name Execution control
2.227 + /// The only way to execute the algorithm is to call the run()
2.228 + /// function.
2.229 +
2.230 + /// @{
2.231 +
2.232 /// \brief Runs the algorithm.
2.233 ///
2.234 /// Runs the algorithm.
2.235 ///
2.236 /// \return \c true if a feasible flow can be found.
2.237 bool run() {
2.238 - init() && start();
2.239 + return init() && start();
2.240 }
2.241
2.242 + /// @}
2.243 +
2.244 + /// \name Query Functions
2.245 + /// The result of the algorithm can be obtained using these
2.246 + /// functions.
2.247 + /// \n run() must be called before using them.
2.248 +
2.249 + /// @{
2.250 +
2.251 /// \brief Returns a const reference to the edge map storing the
2.252 /// found flow.
2.253 ///
2.254 @@ -334,7 +385,7 @@
2.255 ///
2.256 /// \pre \ref run() must be called before using this function.
2.257 const FlowMap& flowMap() const {
2.258 - return _flow;
2.259 + return *_flow;
2.260 }
2.261
2.262 /// \brief Returns a const reference to the node map storing the
2.263 @@ -345,7 +396,25 @@
2.264 ///
2.265 /// \pre \ref run() must be called before using this function.
2.266 const PotentialMap& potentialMap() const {
2.267 - return _potential;
2.268 + return *_potential;
2.269 + }
2.270 +
2.271 + /// \brief Returns the flow on the edge.
2.272 + ///
2.273 + /// Returns the flow on the edge.
2.274 + ///
2.275 + /// \pre \ref run() must be called before using this function.
2.276 + Capacity flow(const Edge& edge) const {
2.277 + return (*_flow)[edge];
2.278 + }
2.279 +
2.280 + /// \brief Returns the potential of the node.
2.281 + ///
2.282 + /// Returns the potential of the node.
2.283 + ///
2.284 + /// \pre \ref run() must be called before using this function.
2.285 + Cost potential(const Node& node) const {
2.286 + return (*_potential)[node];
2.287 }
2.288
2.289 /// \brief Returns the total cost of the found flow.
2.290 @@ -357,16 +426,31 @@
2.291 Cost totalCost() const {
2.292 Cost c = 0;
2.293 for (EdgeIt e(_graph); e != INVALID; ++e)
2.294 - c += _flow[e] * _orig_cost[e];
2.295 + c += (*_flow)[e] * _orig_cost[e];
2.296 return c;
2.297 }
2.298
2.299 + /// @}
2.300 +
2.301 private:
2.302
2.303 /// Initializes the algorithm.
2.304 bool init() {
2.305 if (!_valid_supply) return false;
2.306
2.307 + // Initializing flow and potential maps
2.308 + if (!_flow) {
2.309 + _flow = new FlowMap(_graph);
2.310 + _local_flow = true;
2.311 + }
2.312 + if (!_potential) {
2.313 + _potential = new PotentialMap(_graph);
2.314 + _local_potential = true;
2.315 + }
2.316 +
2.317 + _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
2.318 + _res_graph = new ResGraph(_graph, _capacity, *_flow);
2.319 +
2.320 // Initializing the scaled cost map and the epsilon parameter
2.321 Cost max_cost = 0;
2.322 int node_num = countNodes(_graph);
2.323 @@ -379,9 +463,9 @@
2.324 // Finding a feasible flow using Circulation
2.325 Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
2.326 SupplyMap >
2.327 - circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
2.328 + circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
2.329 _supply );
2.330 - return circulation.flowMap(_flow).run();
2.331 + return circulation.flowMap(*_flow).run();
2.332 }
2.333
2.334
2.335 @@ -397,9 +481,9 @@
2.336 // Performing price refinement heuristic using Bellman-Ford
2.337 // algorithm
2.338 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
2.339 - typedef ShiftMap<ResidualCostMap> ShiftCostMap;
2.340 + typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
2.341 ShiftCostMap shift_cost(_res_cost, _epsilon);
2.342 - BellmanFord<ResGraph, ShiftCostMap> bf(_res_graph, shift_cost);
2.343 + BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
2.344 bf.init(0);
2.345 bool done = false;
2.346 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
2.347 @@ -407,7 +491,7 @@
2.348 done = bf.processNextWeakRound();
2.349 if (done) {
2.350 for (NodeIt n(_graph); n != INVALID; ++n)
2.351 - _potential[n] = bf.dist(n);
2.352 + (*_potential)[n] = bf.dist(n);
2.353 continue;
2.354 }
2.355 }
2.356 @@ -415,16 +499,16 @@
2.357 // Saturating edges not satisfying the optimality condition
2.358 Capacity delta;
2.359 for (EdgeIt e(_graph); e != INVALID; ++e) {
2.360 - if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
2.361 - delta = _capacity[e] - _flow[e];
2.362 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
2.363 + delta = _capacity[e] - (*_flow)[e];
2.364 _excess[_graph.source(e)] -= delta;
2.365 _excess[_graph.target(e)] += delta;
2.366 - _flow[e] = _capacity[e];
2.367 + (*_flow)[e] = _capacity[e];
2.368 }
2.369 - if (_flow[e] > 0 && -_red_cost[e] < 0) {
2.370 - _excess[_graph.target(e)] -= _flow[e];
2.371 - _excess[_graph.source(e)] += _flow[e];
2.372 - _flow[e] = 0;
2.373 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
2.374 + _excess[_graph.target(e)] -= (*_flow)[e];
2.375 + _excess[_graph.source(e)] += (*_flow)[e];
2.376 + (*_flow)[e] = 0;
2.377 }
2.378 }
2.379
2.380 @@ -440,26 +524,26 @@
2.381 // Performing push operations if there are admissible edges
2.382 if (_excess[n] > 0) {
2.383 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
2.384 - if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
2.385 - delta = _capacity[e] - _flow[e] <= _excess[n] ?
2.386 - _capacity[e] - _flow[e] : _excess[n];
2.387 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
2.388 + delta = _capacity[e] - (*_flow)[e] <= _excess[n] ?
2.389 + _capacity[e] - (*_flow)[e] : _excess[n];
2.390 t = _graph.target(e);
2.391
2.392 // Push-look-ahead heuristic
2.393 Capacity ahead = -_excess[t];
2.394 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
2.395 - if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
2.396 - ahead += _capacity[oe] - _flow[oe];
2.397 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
2.398 + ahead += _capacity[oe] - (*_flow)[oe];
2.399 }
2.400 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
2.401 - if (_flow[ie] > 0 && -_red_cost[ie] < 0)
2.402 - ahead += _flow[ie];
2.403 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
2.404 + ahead += (*_flow)[ie];
2.405 }
2.406 if (ahead < 0) ahead = 0;
2.407
2.408 // Pushing flow along the edge
2.409 if (ahead < delta) {
2.410 - _flow[e] += ahead;
2.411 + (*_flow)[e] += ahead;
2.412 _excess[n] -= ahead;
2.413 _excess[t] += ahead;
2.414 active_nodes.push_front(t);
2.415 @@ -467,7 +551,7 @@
2.416 relabel_enabled = false;
2.417 break;
2.418 } else {
2.419 - _flow[e] += delta;
2.420 + (*_flow)[e] += delta;
2.421 _excess[n] -= delta;
2.422 _excess[t] += delta;
2.423 if (_excess[t] > 0 && _excess[t] <= delta)
2.424 @@ -481,25 +565,25 @@
2.425
2.426 if (_excess[n] > 0) {
2.427 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
2.428 - if (_flow[e] > 0 && -_red_cost[e] < 0) {
2.429 - delta = _flow[e] <= _excess[n] ? _flow[e] : _excess[n];
2.430 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
2.431 + delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];
2.432 t = _graph.source(e);
2.433
2.434 // Push-look-ahead heuristic
2.435 Capacity ahead = -_excess[t];
2.436 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
2.437 - if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
2.438 - ahead += _capacity[oe] - _flow[oe];
2.439 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
2.440 + ahead += _capacity[oe] - (*_flow)[oe];
2.441 }
2.442 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
2.443 - if (_flow[ie] > 0 && -_red_cost[ie] < 0)
2.444 - ahead += _flow[ie];
2.445 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
2.446 + ahead += (*_flow)[ie];
2.447 }
2.448 if (ahead < 0) ahead = 0;
2.449
2.450 // Pushing flow along the edge
2.451 if (ahead < delta) {
2.452 - _flow[e] -= ahead;
2.453 + (*_flow)[e] -= ahead;
2.454 _excess[n] -= ahead;
2.455 _excess[t] += ahead;
2.456 active_nodes.push_front(t);
2.457 @@ -507,7 +591,7 @@
2.458 relabel_enabled = false;
2.459 break;
2.460 } else {
2.461 - _flow[e] -= delta;
2.462 + (*_flow)[e] -= delta;
2.463 _excess[n] -= delta;
2.464 _excess[t] += delta;
2.465 if (_excess[t] > 0 && _excess[t] <= delta)
2.466 @@ -523,15 +607,15 @@
2.467 // Performing relabel operation if the node is still active
2.468 LCost min_red_cost = std::numeric_limits<LCost>::max();
2.469 for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
2.470 - if ( _capacity[oe] - _flow[oe] > 0 &&
2.471 - _red_cost[oe] < min_red_cost )
2.472 - min_red_cost = _red_cost[oe];
2.473 + if ( _capacity[oe] - (*_flow)[oe] > 0 &&
2.474 + (*_red_cost)[oe] < min_red_cost )
2.475 + min_red_cost = (*_red_cost)[oe];
2.476 }
2.477 for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
2.478 - if (_flow[ie] > 0 && -_red_cost[ie] < min_red_cost)
2.479 - min_red_cost = -_red_cost[ie];
2.480 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
2.481 + min_red_cost = -(*_red_cost)[ie];
2.482 }
2.483 - _potential[n] -= min_red_cost + _epsilon;
2.484 + (*_potential)[n] -= min_red_cost + _epsilon;
2.485 hyper[n] = false;
2.486 }
2.487
2.488 @@ -544,10 +628,18 @@
2.489 }
2.490 }
2.491
2.492 + // Computing node potentials for the original costs
2.493 + ResidualCostMap<CostMap> res_cost(_orig_cost);
2.494 + BellmanFord< ResGraph, ResidualCostMap<CostMap> >
2.495 + bf(*_res_graph, res_cost);
2.496 + bf.init(0); bf.start();
2.497 + for (NodeIt n(_graph); n != INVALID; ++n)
2.498 + (*_potential)[n] = bf.dist(n);
2.499 +
2.500 // Handling non-zero lower bounds
2.501 if (_lower) {
2.502 for (EdgeIt e(_graph); e != INVALID; ++e)
2.503 - _flow[e] += (*_lower)[e];
2.504 + (*_flow)[e] += (*_lower)[e];
2.505 }
2.506 return true;
2.507 }
3.1 --- a/lemon/cycle_canceling.h Wed Feb 27 11:39:03 2008 +0000
3.2 +++ b/lemon/cycle_canceling.h Thu Feb 28 02:54:27 2008 +0000
3.3 @@ -52,11 +52,8 @@
3.4 /// \warning
3.5 /// - Edge capacities and costs should be \e non-negative \e integers.
3.6 /// - Supply values should be \e signed \e integers.
3.7 - /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
3.8 - /// - \c CapacityMap::Value and \c SupplyMap::Value must be
3.9 - /// convertible to each other.
3.10 - /// - All value types must be convertible to \c CostMap::Value, which
3.11 - /// must be signed type.
3.12 + /// - The value types of the maps should be convertible to each other.
3.13 + /// - \c CostMap::Value must be signed type.
3.14 ///
3.15 /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
3.16 /// used for negative cycle detection with limited iteration number.
3.17 @@ -94,6 +91,8 @@
3.18
3.19 /// The type of the flow map.
3.20 typedef typename Graph::template EdgeMap<Capacity> FlowMap;
3.21 + /// The type of the potential map.
3.22 + typedef typename Graph::template NodeMap<Cost> PotentialMap;
3.23
3.24 private:
3.25
3.26 @@ -143,18 +142,22 @@
3.27 bool _valid_supply;
3.28
3.29 // Edge map of the current flow
3.30 - FlowMap _flow;
3.31 + FlowMap *_flow;
3.32 + bool _local_flow;
3.33 + // Node map of the current potentials
3.34 + PotentialMap *_potential;
3.35 + bool _local_potential;
3.36
3.37 // The residual graph
3.38 - ResGraph _res_graph;
3.39 + ResGraph *_res_graph;
3.40 // The residual cost map
3.41 ResidualCostMap _res_cost;
3.42
3.43 public:
3.44
3.45 - /// \brief General constructor of the class (with lower bounds).
3.46 + /// \brief General constructor (with lower bounds).
3.47 ///
3.48 - /// General constructor of the class (with lower bounds).
3.49 + /// General constructor (with lower bounds).
3.50 ///
3.51 /// \param graph The directed graph the algorithm runs on.
3.52 /// \param lower The lower bounds of the edges.
3.53 @@ -167,8 +170,8 @@
3.54 const CostMap &cost,
3.55 const SupplyMap &supply ) :
3.56 _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
3.57 - _supply(graph), _flow(graph, 0),
3.58 - _res_graph(graph, _capacity, _flow), _res_cost(_cost)
3.59 + _supply(graph), _flow(0), _local_flow(false),
3.60 + _potential(0), _local_potential(false), _res_cost(_cost)
3.61 {
3.62 // Removing non-zero lower bounds
3.63 _capacity = subMap(capacity, lower);
3.64 @@ -184,9 +187,9 @@
3.65 _valid_supply = sum == 0;
3.66 }
3.67
3.68 - /// \brief General constructor of the class (without lower bounds).
3.69 + /// \brief General constructor (without lower bounds).
3.70 ///
3.71 - /// General constructor of the class (without lower bounds).
3.72 + /// General constructor (without lower bounds).
3.73 ///
3.74 /// \param graph The directed graph the algorithm runs on.
3.75 /// \param capacity The capacities (upper bounds) of the edges.
3.76 @@ -197,8 +200,8 @@
3.77 const CostMap &cost,
3.78 const SupplyMap &supply ) :
3.79 _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
3.80 - _supply(supply), _flow(graph, 0),
3.81 - _res_graph(graph, _capacity, _flow), _res_cost(_cost)
3.82 + _supply(supply), _flow(0), _local_flow(false),
3.83 + _potential(0), _local_potential(false), _res_cost(_cost)
3.84 {
3.85 // Checking the sum of supply values
3.86 Supply sum = 0;
3.87 @@ -206,9 +209,9 @@
3.88 _valid_supply = sum == 0;
3.89 }
3.90
3.91 - /// \brief Simple constructor of the class (with lower bounds).
3.92 + /// \brief Simple constructor (with lower bounds).
3.93 ///
3.94 - /// Simple constructor of the class (with lower bounds).
3.95 + /// Simple constructor (with lower bounds).
3.96 ///
3.97 /// \param graph The directed graph the algorithm runs on.
3.98 /// \param lower The lower bounds of the edges.
3.99 @@ -225,8 +228,8 @@
3.100 Node s, Node t,
3.101 Supply flow_value ) :
3.102 _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
3.103 - _supply(graph), _flow(graph, 0),
3.104 - _res_graph(graph, _capacity, _flow), _res_cost(_cost)
3.105 + _supply(graph), _flow(0), _local_flow(false),
3.106 + _potential(0), _local_potential(false), _res_cost(_cost)
3.107 {
3.108 // Removing non-zero lower bounds
3.109 _capacity = subMap(capacity, lower);
3.110 @@ -243,9 +246,9 @@
3.111 _valid_supply = true;
3.112 }
3.113
3.114 - /// \brief Simple constructor of the class (without lower bounds).
3.115 + /// \brief Simple constructor (without lower bounds).
3.116 ///
3.117 - /// Simple constructor of the class (without lower bounds).
3.118 + /// Simple constructor (without lower bounds).
3.119 ///
3.120 /// \param graph The directed graph the algorithm runs on.
3.121 /// \param capacity The capacities (upper bounds) of the edges.
3.122 @@ -260,14 +263,55 @@
3.123 Node s, Node t,
3.124 Supply flow_value ) :
3.125 _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
3.126 - _supply(graph, 0), _flow(graph, 0),
3.127 - _res_graph(graph, _capacity, _flow), _res_cost(_cost)
3.128 + _supply(graph, 0), _flow(0), _local_flow(false),
3.129 + _potential(0), _local_potential(false), _res_cost(_cost)
3.130 {
3.131 _supply[s] = flow_value;
3.132 _supply[t] = -flow_value;
3.133 _valid_supply = true;
3.134 }
3.135
3.136 + /// Destructor.
3.137 + ~CycleCanceling() {
3.138 + if (_local_flow) delete _flow;
3.139 + if (_local_potential) delete _potential;
3.140 + delete _res_graph;
3.141 + }
3.142 +
3.143 + /// \brief Sets the flow map.
3.144 + ///
3.145 + /// Sets the flow map.
3.146 + ///
3.147 + /// \return \c (*this)
3.148 + CycleCanceling& flowMap(FlowMap &map) {
3.149 + if (_local_flow) {
3.150 + delete _flow;
3.151 + _local_flow = false;
3.152 + }
3.153 + _flow = ↦
3.154 + return *this;
3.155 + }
3.156 +
3.157 + /// \brief Sets the potential map.
3.158 + ///
3.159 + /// Sets the potential map.
3.160 + ///
3.161 + /// \return \c (*this)
3.162 + CycleCanceling& potentialMap(PotentialMap &map) {
3.163 + if (_local_potential) {
3.164 + delete _potential;
3.165 + _local_potential = false;
3.166 + }
3.167 + _potential = ↦
3.168 + return *this;
3.169 + }
3.170 +
3.171 + /// \name Execution control
3.172 + /// The only way to execute the algorithm is to call the run()
3.173 + /// function.
3.174 +
3.175 + /// @{
3.176 +
3.177 /// \brief Runs the algorithm.
3.178 ///
3.179 /// Runs the algorithm.
3.180 @@ -281,6 +325,15 @@
3.181 return init() && start(min_mean_cc);
3.182 }
3.183
3.184 + /// @}
3.185 +
3.186 + /// \name Query Functions
3.187 + /// The result of the algorithm can be obtained using these
3.188 + /// functions.
3.189 + /// \n run() must be called before using them.
3.190 +
3.191 + /// @{
3.192 +
3.193 /// \brief Returns a const reference to the edge map storing the
3.194 /// found flow.
3.195 ///
3.196 @@ -288,7 +341,36 @@
3.197 ///
3.198 /// \pre \ref run() must be called before using this function.
3.199 const FlowMap& flowMap() const {
3.200 - return _flow;
3.201 + return *_flow;
3.202 + }
3.203 +
3.204 + /// \brief Returns a const reference to the node map storing the
3.205 + /// found potentials (the dual solution).
3.206 + ///
3.207 + /// Returns a const reference to the node map storing the found
3.208 + /// potentials (the dual solution).
3.209 + ///
3.210 + /// \pre \ref run() must be called before using this function.
3.211 + const PotentialMap& potentialMap() const {
3.212 + return *_potential;
3.213 + }
3.214 +
3.215 + /// \brief Returns the flow on the edge.
3.216 + ///
3.217 + /// Returns the flow on the edge.
3.218 + ///
3.219 + /// \pre \ref run() must be called before using this function.
3.220 + Capacity flow(const Edge& edge) const {
3.221 + return (*_flow)[edge];
3.222 + }
3.223 +
3.224 + /// \brief Returns the potential of the node.
3.225 + ///
3.226 + /// Returns the potential of the node.
3.227 + ///
3.228 + /// \pre \ref run() must be called before using this function.
3.229 + Cost potential(const Node& node) const {
3.230 + return (*_potential)[node];
3.231 }
3.232
3.233 /// \brief Returns the total cost of the found flow.
3.234 @@ -300,29 +382,50 @@
3.235 Cost totalCost() const {
3.236 Cost c = 0;
3.237 for (EdgeIt e(_graph); e != INVALID; ++e)
3.238 - c += _flow[e] * _cost[e];
3.239 + c += (*_flow)[e] * _cost[e];
3.240 return c;
3.241 }
3.242
3.243 + /// @}
3.244 +
3.245 private:
3.246
3.247 /// Initializes the algorithm.
3.248 bool init() {
3.249 if (!_valid_supply) return false;
3.250
3.251 + // Initializing flow and potential maps
3.252 + if (!_flow) {
3.253 + _flow = new FlowMap(_graph);
3.254 + _local_flow = true;
3.255 + }
3.256 + if (!_potential) {
3.257 + _potential = new PotentialMap(_graph);
3.258 + _local_potential = true;
3.259 + }
3.260 +
3.261 + _res_graph = new ResGraph(_graph, _capacity, *_flow);
3.262 +
3.263 // Finding a feasible flow using Circulation
3.264 Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
3.265 SupplyMap >
3.266 - circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
3.267 + circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
3.268 _supply );
3.269 - return circulation.flowMap(_flow).run();
3.270 + return circulation.flowMap(*_flow).run();
3.271 }
3.272
3.273 bool start(bool min_mean_cc) {
3.274 if (min_mean_cc)
3.275 - return startMinMean();
3.276 + startMinMean();
3.277 else
3.278 - return start();
3.279 + start();
3.280 +
3.281 + // Handling non-zero lower bounds
3.282 + if (_lower) {
3.283 + for (EdgeIt e(_graph); e != INVALID; ++e)
3.284 + (*_flow)[e] += (*_lower)[e];
3.285 + }
3.286 + return true;
3.287 }
3.288
3.289 /// \brief Executes the algorithm using \ref BellmanFord.
3.290 @@ -330,16 +433,16 @@
3.291 /// Executes the algorithm using the \ref BellmanFord
3.292 /// "Bellman-Ford" algorithm for negative cycle detection with
3.293 /// successively larger limit for the number of iterations.
3.294 - bool start() {
3.295 - typename BellmanFord<ResGraph, ResidualCostMap>::PredMap pred(_res_graph);
3.296 - typename ResGraph::template NodeMap<int> visited(_res_graph);
3.297 + void start() {
3.298 + typename BellmanFord<ResGraph, ResidualCostMap>::PredMap pred(*_res_graph);
3.299 + typename ResGraph::template NodeMap<int> visited(*_res_graph);
3.300 std::vector<ResEdge> cycle;
3.301 int node_num = countNodes(_graph);
3.302
3.303 int length_bound = BF_FIRST_LIMIT;
3.304 bool optimal = false;
3.305 while (!optimal) {
3.306 - BellmanFord<ResGraph, ResidualCostMap> bf(_res_graph, _res_cost);
3.307 + BellmanFord<ResGraph, ResidualCostMap> bf(*_res_graph, _res_cost);
3.308 bf.predMap(pred);
3.309 bf.init(0);
3.310 int iter_num = 0;
3.311 @@ -356,22 +459,26 @@
3.312 }
3.313 }
3.314 if (real_iter_num < curr_iter_num) {
3.315 + // Optimal flow is found
3.316 optimal = true;
3.317 + // Setting node potentials
3.318 + for (NodeIt n(_graph); n != INVALID; ++n)
3.319 + (*_potential)[n] = bf.dist(n);
3.320 break;
3.321 } else {
3.322 // Searching for node disjoint negative cycles
3.323 - for (ResNodeIt n(_res_graph); n != INVALID; ++n)
3.324 + for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
3.325 visited[n] = 0;
3.326 int id = 0;
3.327 - for (ResNodeIt n(_res_graph); n != INVALID; ++n) {
3.328 + for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
3.329 if (visited[n] > 0) continue;
3.330 visited[n] = ++id;
3.331 ResNode u = pred[n] == INVALID ?
3.332 - INVALID : _res_graph.source(pred[n]);
3.333 + INVALID : _res_graph->source(pred[n]);
3.334 while (u != INVALID && visited[u] == 0) {
3.335 visited[u] = id;
3.336 u = pred[u] == INVALID ?
3.337 - INVALID : _res_graph.source(pred[u]);
3.338 + INVALID : _res_graph->source(pred[u]);
3.339 }
3.340 if (u != INVALID && visited[u] == id) {
3.341 // Finding the negative cycle
3.342 @@ -379,16 +486,16 @@
3.343 cycle.clear();
3.344 ResEdge e = pred[u];
3.345 cycle.push_back(e);
3.346 - Capacity d = _res_graph.rescap(e);
3.347 - while (_res_graph.source(e) != u) {
3.348 - cycle.push_back(e = pred[_res_graph.source(e)]);
3.349 - if (_res_graph.rescap(e) < d)
3.350 - d = _res_graph.rescap(e);
3.351 + Capacity d = _res_graph->rescap(e);
3.352 + while (_res_graph->source(e) != u) {
3.353 + cycle.push_back(e = pred[_res_graph->source(e)]);
3.354 + if (_res_graph->rescap(e) < d)
3.355 + d = _res_graph->rescap(e);
3.356 }
3.357
3.358 // Augmenting along the cycle
3.359 for (int i = 0; i < int(cycle.size()); ++i)
3.360 - _res_graph.augment(cycle[i], d);
3.361 + _res_graph->augment(cycle[i], d);
3.362 }
3.363 }
3.364 }
3.365 @@ -397,22 +504,15 @@
3.366 length_bound = int(length_bound * BF_ALPHA);
3.367 }
3.368 }
3.369 -
3.370 - // Handling non-zero lower bounds
3.371 - if (_lower) {
3.372 - for (EdgeIt e(_graph); e != INVALID; ++e)
3.373 - _flow[e] += (*_lower)[e];
3.374 - }
3.375 - return true;
3.376 }
3.377
3.378 /// \brief Executes the algorithm using \ref MinMeanCycle.
3.379 ///
3.380 /// Executes the algorithm using \ref MinMeanCycle for negative
3.381 /// cycle detection.
3.382 - bool startMinMean() {
3.383 + void startMinMean() {
3.384 typedef Path<ResGraph> ResPath;
3.385 - MinMeanCycle<ResGraph, ResidualCostMap> mmc(_res_graph, _res_cost);
3.386 + MinMeanCycle<ResGraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
3.387 ResPath cycle;
3.388
3.389 mmc.cyclePath(cycle).init();
3.390 @@ -425,13 +525,13 @@
3.391 // along the cycle
3.392 Capacity delta = 0;
3.393 for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) {
3.394 - if (delta == 0 || _res_graph.rescap(e) < delta)
3.395 - delta = _res_graph.rescap(e);
3.396 + if (delta == 0 || _res_graph->rescap(e) < delta)
3.397 + delta = _res_graph->rescap(e);
3.398 }
3.399
3.400 // Augmenting along the cycle
3.401 for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e)
3.402 - _res_graph.augment(e, delta);
3.403 + _res_graph->augment(e, delta);
3.404
3.405 // Finding the minimum cycle mean for the modified residual
3.406 // graph
3.407 @@ -440,12 +540,11 @@
3.408 }
3.409 }
3.410
3.411 - // Handling non-zero lower bounds
3.412 - if (_lower) {
3.413 - for (EdgeIt e(_graph); e != INVALID; ++e)
3.414 - _flow[e] += (*_lower)[e];
3.415 - }
3.416 - return true;
3.417 + // Computing node potentials
3.418 + BellmanFord<ResGraph, ResidualCostMap> bf(*_res_graph, _res_cost);
3.419 + bf.init(0); bf.start();
3.420 + for (NodeIt n(_graph); n != INVALID; ++n)
3.421 + (*_potential)[n] = bf.dist(n);
3.422 }
3.423
3.424 }; //class CycleCanceling
4.1 --- a/lemon/min_cost_flow.h Wed Feb 27 11:39:03 2008 +0000
4.2 +++ b/lemon/min_cost_flow.h Thu Feb 28 02:54:27 2008 +0000
4.3 @@ -43,9 +43,7 @@
4.4 /// \ref NetworkSimplex.
4.5 ///
4.6 /// There are four implementations for the minimum cost flow problem,
4.7 - /// which can be used exactly the same way except for the fact that
4.8 - /// \ref CycleCanceling does not provide dual solution (node
4.9 - /// potentials) since it is a pure primal method.
4.10 + /// which can be used exactly the same way.
4.11 /// - \ref CapacityScaling The capacity scaling algorithm.
4.12 /// - \ref CostScaling The cost scaling algorithm.
4.13 /// - \ref CycleCanceling A cycle-canceling algorithm.
4.14 @@ -60,20 +58,16 @@
4.15 /// \warning
4.16 /// - Edge capacities and costs should be \e non-negative \e integers.
4.17 /// - Supply values should be \e signed \e integers.
4.18 - /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
4.19 - /// - \c CapacityMap::Value and \c SupplyMap::Value must be
4.20 - /// convertible to each other.
4.21 - /// - All value types must be convertible to \c CostMap::Value, which
4.22 - /// must be signed type.
4.23 + /// - The value types of the maps should be convertible to each other.
4.24 + /// - \c CostMap::Value must be signed type.
4.25 ///
4.26 /// \author Peter Kovacs
4.27
4.28 template < typename Graph,
4.29 typename LowerMap = typename Graph::template EdgeMap<int>,
4.30 - typename CapacityMap = LowerMap,
4.31 + typename CapacityMap = typename Graph::template EdgeMap<int>,
4.32 typename CostMap = typename Graph::template EdgeMap<int>,
4.33 - typename SupplyMap = typename Graph::template NodeMap
4.34 - <typename CapacityMap::Value> >
4.35 + typename SupplyMap = typename Graph::template NodeMap<int> >
4.36 class MinCostFlow :
4.37 public NetworkSimplex< Graph, LowerMap, CapacityMap,
4.38 CostMap, SupplyMap >
4.39 @@ -85,7 +79,7 @@
4.40
4.41 public:
4.42
4.43 - /// General constructor of the class (with lower bounds).
4.44 + /// General constructor (with lower bounds).
4.45 MinCostFlow( const Graph &graph,
4.46 const LowerMap &lower,
4.47 const CapacityMap &capacity,
4.48 @@ -100,7 +94,7 @@
4.49 const SupplyMap &supply ) :
4.50 MinCostFlowImpl(graph, capacity, cost, supply) {}
4.51
4.52 - /// Simple constructor of the class (with lower bounds).
4.53 + /// Simple constructor (with lower bounds).
4.54 MinCostFlow( const Graph &graph,
4.55 const LowerMap &lower,
4.56 const CapacityMap &capacity,
4.57 @@ -110,7 +104,7 @@
4.58 MinCostFlowImpl( graph, lower, capacity, cost,
4.59 s, t, flow_value ) {}
4.60
4.61 - /// Simple constructor of the class (without lower bounds).
4.62 + /// Simple constructor (without lower bounds).
4.63 MinCostFlow( const Graph &graph,
4.64 const CapacityMap &capacity,
4.65 const CostMap &cost,
5.1 --- a/lemon/network_simplex.h Wed Feb 27 11:39:03 2008 +0000
5.2 +++ b/lemon/network_simplex.h Thu Feb 28 02:54:27 2008 +0000
5.3 @@ -52,11 +52,8 @@
5.4 /// \warning
5.5 /// - Edge capacities and costs should be \e non-negative \e integers.
5.6 /// - Supply values should be \e signed \e integers.
5.7 - /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
5.8 - /// - \c CapacityMap::Value and \c SupplyMap::Value must be
5.9 - /// convertible to each other.
5.10 - /// - All value types must be convertible to \c CostMap::Value, which
5.11 - /// must be signed type.
5.12 + /// - The value types of the maps should be convertible to each other.
5.13 + /// - \c CostMap::Value must be signed type.
5.14 ///
5.15 /// \note \ref NetworkSimplex provides six different pivot rule
5.16 /// implementations that significantly affect the efficiency of the
5.17 @@ -469,8 +466,10 @@
5.18 ReducedCostMap _red_cost;
5.19
5.20 // Members for handling the original graph
5.21 - FlowMap _flow_result;
5.22 - PotentialMap _potential_result;
5.23 + FlowMap *_flow_result;
5.24 + PotentialMap *_potential_result;
5.25 + bool _local_flow;
5.26 + bool _local_potential;
5.27 NodeRefMap _node_ref;
5.28 EdgeRefMap _edge_ref;
5.29
5.30 @@ -487,9 +486,9 @@
5.31
5.32 public :
5.33
5.34 - /// \brief General constructor of the class (with lower bounds).
5.35 + /// \brief General constructor (with lower bounds).
5.36 ///
5.37 - /// General constructor of the class (with lower bounds).
5.38 + /// General constructor (with lower bounds).
5.39 ///
5.40 /// \param graph The directed graph the algorithm runs on.
5.41 /// \param lower The lower bounds of the edges.
5.42 @@ -506,7 +505,8 @@
5.43 _potential(_graph), _depth(_graph), _parent(_graph),
5.44 _pred_edge(_graph), _thread(_graph), _forward(_graph),
5.45 _state(_graph), _red_cost(_graph, _cost, _potential),
5.46 - _flow_result(graph), _potential_result(graph),
5.47 + _flow_result(0), _potential_result(0),
5.48 + _local_flow(false), _local_potential(false),
5.49 _node_ref(graph), _edge_ref(graph)
5.50 {
5.51 // Checking the sum of supply values
5.52 @@ -538,9 +538,9 @@
5.53 }
5.54 }
5.55
5.56 - /// \brief General constructor of the class (without lower bounds).
5.57 + /// \brief General constructor (without lower bounds).
5.58 ///
5.59 - /// General constructor of the class (without lower bounds).
5.60 + /// General constructor (without lower bounds).
5.61 ///
5.62 /// \param graph The directed graph the algorithm runs on.
5.63 /// \param capacity The capacities (upper bounds) of the edges.
5.64 @@ -555,7 +555,8 @@
5.65 _potential(_graph), _depth(_graph), _parent(_graph),
5.66 _pred_edge(_graph), _thread(_graph), _forward(_graph),
5.67 _state(_graph), _red_cost(_graph, _cost, _potential),
5.68 - _flow_result(graph), _potential_result(graph),
5.69 + _flow_result(0), _potential_result(0),
5.70 + _local_flow(false), _local_potential(false),
5.71 _node_ref(graph), _edge_ref(graph)
5.72 {
5.73 // Checking the sum of supply values
5.74 @@ -574,9 +575,9 @@
5.75 .run();
5.76 }
5.77
5.78 - /// \brief Simple constructor of the class (with lower bounds).
5.79 + /// \brief Simple constructor (with lower bounds).
5.80 ///
5.81 - /// Simple constructor of the class (with lower bounds).
5.82 + /// Simple constructor (with lower bounds).
5.83 ///
5.84 /// \param graph The directed graph the algorithm runs on.
5.85 /// \param lower The lower bounds of the edges.
5.86 @@ -598,7 +599,8 @@
5.87 _potential(_graph), _depth(_graph), _parent(_graph),
5.88 _pred_edge(_graph), _thread(_graph), _forward(_graph),
5.89 _state(_graph), _red_cost(_graph, _cost, _potential),
5.90 - _flow_result(graph), _potential_result(graph),
5.91 + _flow_result(0), _potential_result(0),
5.92 + _local_flow(false), _local_potential(false),
5.93 _node_ref(graph), _edge_ref(graph)
5.94 {
5.95 // Copying _graph_ref to graph
5.96 @@ -625,9 +627,9 @@
5.97 _valid_supply = true;
5.98 }
5.99
5.100 - /// \brief Simple constructor of the class (without lower bounds).
5.101 + /// \brief Simple constructor (without lower bounds).
5.102 ///
5.103 - /// Simple constructor of the class (without lower bounds).
5.104 + /// Simple constructor (without lower bounds).
5.105 ///
5.106 /// \param graph The directed graph the algorithm runs on.
5.107 /// \param capacity The capacities (upper bounds) of the edges.
5.108 @@ -647,7 +649,8 @@
5.109 _potential(_graph), _depth(_graph), _parent(_graph),
5.110 _pred_edge(_graph), _thread(_graph), _forward(_graph),
5.111 _state(_graph), _red_cost(_graph, _cost, _potential),
5.112 - _flow_result(graph), _potential_result(graph),
5.113 + _flow_result(0), _potential_result(0),
5.114 + _local_flow(false), _local_potential(false),
5.115 _node_ref(graph), _edge_ref(graph)
5.116 {
5.117 // Copying _graph_ref to graph
5.118 @@ -662,6 +665,46 @@
5.119 _valid_supply = true;
5.120 }
5.121
5.122 + /// Destructor.
5.123 + ~NetworkSimplex() {
5.124 + if (_local_flow) delete _flow_result;
5.125 + if (_local_potential) delete _potential_result;
5.126 + }
5.127 +
5.128 + /// \brief Sets the flow map.
5.129 + ///
5.130 + /// Sets the flow map.
5.131 + ///
5.132 + /// \return \c (*this)
5.133 + NetworkSimplex& flowMap(FlowMap &map) {
5.134 + if (_local_flow) {
5.135 + delete _flow_result;
5.136 + _local_flow = false;
5.137 + }
5.138 + _flow_result = ↦
5.139 + return *this;
5.140 + }
5.141 +
5.142 + /// \brief Sets the potential map.
5.143 + ///
5.144 + /// Sets the potential map.
5.145 + ///
5.146 + /// \return \c (*this)
5.147 + NetworkSimplex& potentialMap(PotentialMap &map) {
5.148 + if (_local_potential) {
5.149 + delete _potential_result;
5.150 + _local_potential = false;
5.151 + }
5.152 + _potential_result = ↦
5.153 + return *this;
5.154 + }
5.155 +
5.156 + /// \name Execution control
5.157 + /// The only way to execute the algorithm is to call the run()
5.158 + /// function.
5.159 +
5.160 + /// @{
5.161 +
5.162 /// \brief Runs the algorithm.
5.163 ///
5.164 /// Runs the algorithm.
5.165 @@ -703,6 +746,15 @@
5.166 return init() && start(pivot_rule);
5.167 }
5.168
5.169 + /// @}
5.170 +
5.171 + /// \name Query Functions
5.172 + /// The result of the algorithm can be obtained using these
5.173 + /// functions.
5.174 + /// \n run() must be called before using them.
5.175 +
5.176 + /// @{
5.177 +
5.178 /// \brief Returns a const reference to the edge map storing the
5.179 /// found flow.
5.180 ///
5.181 @@ -710,7 +762,7 @@
5.182 ///
5.183 /// \pre \ref run() must be called before using this function.
5.184 const FlowMap& flowMap() const {
5.185 - return _flow_result;
5.186 + return *_flow_result;
5.187 }
5.188
5.189 /// \brief Returns a const reference to the node map storing the
5.190 @@ -721,7 +773,25 @@
5.191 ///
5.192 /// \pre \ref run() must be called before using this function.
5.193 const PotentialMap& potentialMap() const {
5.194 - return _potential_result;
5.195 + return *_potential_result;
5.196 + }
5.197 +
5.198 + /// \brief Returns the flow on the edge.
5.199 + ///
5.200 + /// Returns the flow on the edge.
5.201 + ///
5.202 + /// \pre \ref run() must be called before using this function.
5.203 + Capacity flow(const typename Graph::Edge& edge) const {
5.204 + return (*_flow_result)[edge];
5.205 + }
5.206 +
5.207 + /// \brief Returns the potential of the node.
5.208 + ///
5.209 + /// Returns the potential of the node.
5.210 + ///
5.211 + /// \pre \ref run() must be called before using this function.
5.212 + Cost potential(const typename Graph::Node& node) const {
5.213 + return (*_potential_result)[node];
5.214 }
5.215
5.216 /// \brief Returns the total cost of the found flow.
5.217 @@ -733,10 +803,12 @@
5.218 Cost totalCost() const {
5.219 Cost c = 0;
5.220 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
5.221 - c += _flow_result[e] * _cost[_edge_ref[e]];
5.222 + c += (*_flow_result)[e] * _cost[_edge_ref[e]];
5.223 return c;
5.224 }
5.225
5.226 + /// @}
5.227 +
5.228 private:
5.229
5.230 /// \brief Extends the underlaying graph and initializes all the
5.231 @@ -744,6 +816,16 @@
5.232 bool init() {
5.233 if (!_valid_supply) return false;
5.234
5.235 + // Initializing result maps
5.236 + if (!_flow_result) {
5.237 + _flow_result = new FlowMap(_graph_ref);
5.238 + _local_flow = true;
5.239 + }
5.240 + if (!_potential_result) {
5.241 + _potential_result = new PotentialMap(_graph_ref);
5.242 + _local_potential = true;
5.243 + }
5.244 +
5.245 // Initializing state and flow maps
5.246 for (EdgeIt e(_graph); e != INVALID; ++e) {
5.247 _flow[e] = 0;
5.248 @@ -1015,14 +1097,14 @@
5.249 // Copying flow values to _flow_result
5.250 if (_lower) {
5.251 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
5.252 - _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]];
5.253 + (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
5.254 } else {
5.255 for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
5.256 - _flow_result[e] = _flow[_edge_ref[e]];
5.257 + (*_flow_result)[e] = _flow[_edge_ref[e]];
5.258 }
5.259 // Copying potential values to _potential_result
5.260 for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
5.261 - _potential_result[n] = _potential[_node_ref[n]];
5.262 + (*_potential_result)[n] = _potential[_node_ref[n]];
5.263
5.264 return true;
5.265 }