1.1 --- a/lemon/cost_scaling.h Wed Feb 27 11:39:03 2008 +0000
1.2 +++ b/lemon/cost_scaling.h Thu Feb 28 02:54:27 2008 +0000
1.3 @@ -54,11 +54,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 /// \note Edge costs are multiplied with the number of nodes during
1.16 /// the algorithm so overflow problems may arise more easily than with
1.17 @@ -97,7 +94,7 @@
1.18 public:
1.19
1.20 /// The type of the flow map.
1.21 - typedef CapacityEdgeMap FlowMap;
1.22 + typedef typename Graph::template EdgeMap<Capacity> FlowMap;
1.23 /// The type of the potential map.
1.24 typedef typename Graph::template NodeMap<LCost> PotentialMap;
1.25
1.26 @@ -107,20 +104,21 @@
1.27 ///
1.28 /// \ref ResidualCostMap is a map adaptor class for handling
1.29 /// residual edge costs.
1.30 - class ResidualCostMap : public MapBase<ResEdge, LCost>
1.31 + template <typename Map>
1.32 + class ResidualCostMap : public MapBase<ResEdge, typename Map::Value>
1.33 {
1.34 private:
1.35
1.36 - const LargeCostMap &_cost_map;
1.37 + const Map &_cost_map;
1.38
1.39 public:
1.40
1.41 ///\e
1.42 - ResidualCostMap(const LargeCostMap &cost_map) :
1.43 + ResidualCostMap(const Map &cost_map) :
1.44 _cost_map(cost_map) {}
1.45
1.46 ///\e
1.47 - LCost operator[](const ResEdge &e) const {
1.48 + typename Map::Value operator[](const ResEdge &e) const {
1.49 return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.50 }
1.51
1.52 @@ -160,8 +158,8 @@
1.53 static const int ALPHA = 4;
1.54
1.55 // Paramters for heuristics
1.56 - static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
1.57 - static const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.58 + static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
1.59 + static const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.60
1.61 private:
1.62
1.63 @@ -180,16 +178,18 @@
1.64 bool _valid_supply;
1.65
1.66 // Edge map of the current flow
1.67 - FlowMap _flow;
1.68 + FlowMap *_flow;
1.69 + bool _local_flow;
1.70 // Node map of the current potentials
1.71 - PotentialMap _potential;
1.72 + PotentialMap *_potential;
1.73 + bool _local_potential;
1.74
1.75 // The residual graph
1.76 - ResGraph _res_graph;
1.77 + ResGraph *_res_graph;
1.78 // The residual cost map
1.79 - ResidualCostMap _res_cost;
1.80 + ResidualCostMap<LargeCostMap> _res_cost;
1.81 // The reduced cost map
1.82 - ReducedCostMap _red_cost;
1.83 + ReducedCostMap *_red_cost;
1.84 // The excess map
1.85 SupplyNodeMap _excess;
1.86 // The epsilon parameter used for cost scaling
1.87 @@ -197,9 +197,9 @@
1.88
1.89 public:
1.90
1.91 - /// \brief General constructor of the class (with lower bounds).
1.92 + /// \brief General constructor (with lower bounds).
1.93 ///
1.94 - /// General constructor of the class (with lower bounds).
1.95 + /// General constructor (with lower bounds).
1.96 ///
1.97 /// \param graph The directed graph the algorithm runs on.
1.98 /// \param lower The lower bounds of the edges.
1.99 @@ -212,9 +212,9 @@
1.100 const CostMap &cost,
1.101 const SupplyMap &supply ) :
1.102 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
1.103 - _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
1.104 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
1.105 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
1.106 + _cost(graph), _supply(graph), _flow(0), _local_flow(false),
1.107 + _potential(0), _local_potential(false), _res_cost(_cost),
1.108 + _excess(graph, 0)
1.109 {
1.110 // Removing non-zero lower bounds
1.111 _capacity = subMap(capacity, lower);
1.112 @@ -231,9 +231,9 @@
1.113 _valid_supply = sum == 0;
1.114 }
1.115
1.116 - /// \brief General constructor of the class (without lower bounds).
1.117 + /// \brief General constructor (without lower bounds).
1.118 ///
1.119 - /// General constructor of the class (without lower bounds).
1.120 + /// General constructor (without lower bounds).
1.121 ///
1.122 /// \param graph The directed graph the algorithm runs on.
1.123 /// \param capacity The capacities (upper bounds) of the edges.
1.124 @@ -244,9 +244,9 @@
1.125 const CostMap &cost,
1.126 const SupplyMap &supply ) :
1.127 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.128 - _cost(graph), _supply(supply), _flow(graph, 0), _potential(graph, 0),
1.129 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
1.130 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
1.131 + _cost(graph), _supply(supply), _flow(0), _local_flow(false),
1.132 + _potential(0), _local_potential(false), _res_cost(_cost),
1.133 + _excess(graph, 0)
1.134 {
1.135 // Checking the sum of supply values
1.136 Supply sum = 0;
1.137 @@ -254,9 +254,9 @@
1.138 _valid_supply = sum == 0;
1.139 }
1.140
1.141 - /// \brief Simple constructor of the class (with lower bounds).
1.142 + /// \brief Simple constructor (with lower bounds).
1.143 ///
1.144 - /// Simple constructor of the class (with lower bounds).
1.145 + /// Simple constructor (with lower bounds).
1.146 ///
1.147 /// \param graph The directed graph the algorithm runs on.
1.148 /// \param lower The lower bounds of the edges.
1.149 @@ -273,9 +273,9 @@
1.150 Node s, Node t,
1.151 Supply flow_value ) :
1.152 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
1.153 - _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
1.154 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
1.155 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
1.156 + _cost(graph), _supply(graph), _flow(0), _local_flow(false),
1.157 + _potential(0), _local_potential(false), _res_cost(_cost),
1.158 + _excess(graph, 0)
1.159 {
1.160 // Removing nonzero lower bounds
1.161 _capacity = subMap(capacity, lower);
1.162 @@ -292,9 +292,9 @@
1.163 _valid_supply = true;
1.164 }
1.165
1.166 - /// \brief Simple constructor of the class (without lower bounds).
1.167 + /// \brief Simple constructor (without lower bounds).
1.168 ///
1.169 - /// Simple constructor of the class (without lower bounds).
1.170 + /// Simple constructor (without lower bounds).
1.171 ///
1.172 /// \param graph The directed graph the algorithm runs on.
1.173 /// \param capacity The capacities (upper bounds) of the edges.
1.174 @@ -309,24 +309,75 @@
1.175 Node s, Node t,
1.176 Supply flow_value ) :
1.177 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.178 - _cost(graph), _supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
1.179 - _res_graph(graph, _capacity, _flow), _res_cost(_cost),
1.180 - _red_cost(graph, _cost, _potential), _excess(graph, 0)
1.181 + _cost(graph), _supply(graph, 0), _flow(0), _local_flow(false),
1.182 + _potential(0), _local_potential(false), _res_cost(_cost),
1.183 + _excess(graph, 0)
1.184 {
1.185 _supply[s] = flow_value;
1.186 _supply[t] = -flow_value;
1.187 _valid_supply = true;
1.188 }
1.189
1.190 + /// Destructor.
1.191 + ~CostScaling() {
1.192 + if (_local_flow) delete _flow;
1.193 + if (_local_potential) delete _potential;
1.194 + delete _res_graph;
1.195 + delete _red_cost;
1.196 + }
1.197 +
1.198 + /// \brief Sets the flow map.
1.199 + ///
1.200 + /// Sets the flow map.
1.201 + ///
1.202 + /// \return \c (*this)
1.203 + CostScaling& flowMap(FlowMap &map) {
1.204 + if (_local_flow) {
1.205 + delete _flow;
1.206 + _local_flow = false;
1.207 + }
1.208 + _flow = ↦
1.209 + return *this;
1.210 + }
1.211 +
1.212 + /// \brief Sets the potential map.
1.213 + ///
1.214 + /// Sets the potential map.
1.215 + ///
1.216 + /// \return \c (*this)
1.217 + CostScaling& potentialMap(PotentialMap &map) {
1.218 + if (_local_potential) {
1.219 + delete _potential;
1.220 + _local_potential = false;
1.221 + }
1.222 + _potential = ↦
1.223 + return *this;
1.224 + }
1.225 +
1.226 + /// \name Execution control
1.227 + /// The only way to execute the algorithm is to call the run()
1.228 + /// function.
1.229 +
1.230 + /// @{
1.231 +
1.232 /// \brief Runs the algorithm.
1.233 ///
1.234 /// Runs the algorithm.
1.235 ///
1.236 /// \return \c true if a feasible flow can be found.
1.237 bool run() {
1.238 - init() && start();
1.239 + return init() && start();
1.240 }
1.241
1.242 + /// @}
1.243 +
1.244 + /// \name Query Functions
1.245 + /// The result of the algorithm can be obtained using these
1.246 + /// functions.
1.247 + /// \n run() must be called before using them.
1.248 +
1.249 + /// @{
1.250 +
1.251 /// \brief Returns a const reference to the edge map storing the
1.252 /// found flow.
1.253 ///
1.254 @@ -334,7 +385,7 @@
1.255 ///
1.256 /// \pre \ref run() must be called before using this function.
1.257 const FlowMap& flowMap() const {
1.258 - return _flow;
1.259 + return *_flow;
1.260 }
1.261
1.262 /// \brief Returns a const reference to the node map storing the
1.263 @@ -345,7 +396,25 @@
1.264 ///
1.265 /// \pre \ref run() must be called before using this function.
1.266 const PotentialMap& potentialMap() const {
1.267 - return _potential;
1.268 + return *_potential;
1.269 + }
1.270 +
1.271 + /// \brief Returns the flow on the edge.
1.272 + ///
1.273 + /// Returns the flow on the edge.
1.274 + ///
1.275 + /// \pre \ref run() must be called before using this function.
1.276 + Capacity flow(const Edge& edge) const {
1.277 + return (*_flow)[edge];
1.278 + }
1.279 +
1.280 + /// \brief Returns the potential of the node.
1.281 + ///
1.282 + /// Returns the potential of the node.
1.283 + ///
1.284 + /// \pre \ref run() must be called before using this function.
1.285 + Cost potential(const Node& node) const {
1.286 + return (*_potential)[node];
1.287 }
1.288
1.289 /// \brief Returns the total cost of the found flow.
1.290 @@ -357,16 +426,31 @@
1.291 Cost totalCost() const {
1.292 Cost c = 0;
1.293 for (EdgeIt e(_graph); e != INVALID; ++e)
1.294 - c += _flow[e] * _orig_cost[e];
1.295 + c += (*_flow)[e] * _orig_cost[e];
1.296 return c;
1.297 }
1.298
1.299 + /// @}
1.300 +
1.301 private:
1.302
1.303 /// Initializes the algorithm.
1.304 bool init() {
1.305 if (!_valid_supply) return false;
1.306
1.307 + // Initializing flow and potential maps
1.308 + if (!_flow) {
1.309 + _flow = new FlowMap(_graph);
1.310 + _local_flow = true;
1.311 + }
1.312 + if (!_potential) {
1.313 + _potential = new PotentialMap(_graph);
1.314 + _local_potential = true;
1.315 + }
1.316 +
1.317 + _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
1.318 + _res_graph = new ResGraph(_graph, _capacity, *_flow);
1.319 +
1.320 // Initializing the scaled cost map and the epsilon parameter
1.321 Cost max_cost = 0;
1.322 int node_num = countNodes(_graph);
1.323 @@ -379,9 +463,9 @@
1.324 // Finding a feasible flow using Circulation
1.325 Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
1.326 SupplyMap >
1.327 - circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
1.328 + circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
1.329 _supply );
1.330 - return circulation.flowMap(_flow).run();
1.331 + return circulation.flowMap(*_flow).run();
1.332 }
1.333
1.334
1.335 @@ -397,9 +481,9 @@
1.336 // Performing price refinement heuristic using Bellman-Ford
1.337 // algorithm
1.338 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.339 - typedef ShiftMap<ResidualCostMap> ShiftCostMap;
1.340 + typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.341 ShiftCostMap shift_cost(_res_cost, _epsilon);
1.342 - BellmanFord<ResGraph, ShiftCostMap> bf(_res_graph, shift_cost);
1.343 + BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.344 bf.init(0);
1.345 bool done = false;
1.346 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.347 @@ -407,7 +491,7 @@
1.348 done = bf.processNextWeakRound();
1.349 if (done) {
1.350 for (NodeIt n(_graph); n != INVALID; ++n)
1.351 - _potential[n] = bf.dist(n);
1.352 + (*_potential)[n] = bf.dist(n);
1.353 continue;
1.354 }
1.355 }
1.356 @@ -415,16 +499,16 @@
1.357 // Saturating edges not satisfying the optimality condition
1.358 Capacity delta;
1.359 for (EdgeIt e(_graph); e != INVALID; ++e) {
1.360 - if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
1.361 - delta = _capacity[e] - _flow[e];
1.362 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.363 + delta = _capacity[e] - (*_flow)[e];
1.364 _excess[_graph.source(e)] -= delta;
1.365 _excess[_graph.target(e)] += delta;
1.366 - _flow[e] = _capacity[e];
1.367 + (*_flow)[e] = _capacity[e];
1.368 }
1.369 - if (_flow[e] > 0 && -_red_cost[e] < 0) {
1.370 - _excess[_graph.target(e)] -= _flow[e];
1.371 - _excess[_graph.source(e)] += _flow[e];
1.372 - _flow[e] = 0;
1.373 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.374 + _excess[_graph.target(e)] -= (*_flow)[e];
1.375 + _excess[_graph.source(e)] += (*_flow)[e];
1.376 + (*_flow)[e] = 0;
1.377 }
1.378 }
1.379
1.380 @@ -440,26 +524,26 @@
1.381 // Performing push operations if there are admissible edges
1.382 if (_excess[n] > 0) {
1.383 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.384 - if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
1.385 - delta = _capacity[e] - _flow[e] <= _excess[n] ?
1.386 - _capacity[e] - _flow[e] : _excess[n];
1.387 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.388 + delta = _capacity[e] - (*_flow)[e] <= _excess[n] ?
1.389 + _capacity[e] - (*_flow)[e] : _excess[n];
1.390 t = _graph.target(e);
1.391
1.392 // Push-look-ahead heuristic
1.393 Capacity ahead = -_excess[t];
1.394 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
1.395 - if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
1.396 - ahead += _capacity[oe] - _flow[oe];
1.397 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.398 + ahead += _capacity[oe] - (*_flow)[oe];
1.399 }
1.400 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
1.401 - if (_flow[ie] > 0 && -_red_cost[ie] < 0)
1.402 - ahead += _flow[ie];
1.403 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.404 + ahead += (*_flow)[ie];
1.405 }
1.406 if (ahead < 0) ahead = 0;
1.407
1.408 // Pushing flow along the edge
1.409 if (ahead < delta) {
1.410 - _flow[e] += ahead;
1.411 + (*_flow)[e] += ahead;
1.412 _excess[n] -= ahead;
1.413 _excess[t] += ahead;
1.414 active_nodes.push_front(t);
1.415 @@ -467,7 +551,7 @@
1.416 relabel_enabled = false;
1.417 break;
1.418 } else {
1.419 - _flow[e] += delta;
1.420 + (*_flow)[e] += delta;
1.421 _excess[n] -= delta;
1.422 _excess[t] += delta;
1.423 if (_excess[t] > 0 && _excess[t] <= delta)
1.424 @@ -481,25 +565,25 @@
1.425
1.426 if (_excess[n] > 0) {
1.427 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.428 - if (_flow[e] > 0 && -_red_cost[e] < 0) {
1.429 - delta = _flow[e] <= _excess[n] ? _flow[e] : _excess[n];
1.430 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.431 + delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];
1.432 t = _graph.source(e);
1.433
1.434 // Push-look-ahead heuristic
1.435 Capacity ahead = -_excess[t];
1.436 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
1.437 - if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
1.438 - ahead += _capacity[oe] - _flow[oe];
1.439 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.440 + ahead += _capacity[oe] - (*_flow)[oe];
1.441 }
1.442 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
1.443 - if (_flow[ie] > 0 && -_red_cost[ie] < 0)
1.444 - ahead += _flow[ie];
1.445 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.446 + ahead += (*_flow)[ie];
1.447 }
1.448 if (ahead < 0) ahead = 0;
1.449
1.450 // Pushing flow along the edge
1.451 if (ahead < delta) {
1.452 - _flow[e] -= ahead;
1.453 + (*_flow)[e] -= ahead;
1.454 _excess[n] -= ahead;
1.455 _excess[t] += ahead;
1.456 active_nodes.push_front(t);
1.457 @@ -507,7 +591,7 @@
1.458 relabel_enabled = false;
1.459 break;
1.460 } else {
1.461 - _flow[e] -= delta;
1.462 + (*_flow)[e] -= delta;
1.463 _excess[n] -= delta;
1.464 _excess[t] += delta;
1.465 if (_excess[t] > 0 && _excess[t] <= delta)
1.466 @@ -523,15 +607,15 @@
1.467 // Performing relabel operation if the node is still active
1.468 LCost min_red_cost = std::numeric_limits<LCost>::max();
1.469 for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
1.470 - if ( _capacity[oe] - _flow[oe] > 0 &&
1.471 - _red_cost[oe] < min_red_cost )
1.472 - min_red_cost = _red_cost[oe];
1.473 + if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.474 + (*_red_cost)[oe] < min_red_cost )
1.475 + min_red_cost = (*_red_cost)[oe];
1.476 }
1.477 for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
1.478 - if (_flow[ie] > 0 && -_red_cost[ie] < min_red_cost)
1.479 - min_red_cost = -_red_cost[ie];
1.480 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.481 + min_red_cost = -(*_red_cost)[ie];
1.482 }
1.483 - _potential[n] -= min_red_cost + _epsilon;
1.484 + (*_potential)[n] -= min_red_cost + _epsilon;
1.485 hyper[n] = false;
1.486 }
1.487
1.488 @@ -544,10 +628,18 @@
1.489 }
1.490 }
1.491
1.492 + // Computing node potentials for the original costs
1.493 + ResidualCostMap<CostMap> res_cost(_orig_cost);
1.494 + BellmanFord< ResGraph, ResidualCostMap<CostMap> >
1.495 + bf(*_res_graph, res_cost);
1.496 + bf.init(0); bf.start();
1.497 + for (NodeIt n(_graph); n != INVALID; ++n)
1.498 + (*_potential)[n] = bf.dist(n);
1.499 +
1.500 // Handling non-zero lower bounds
1.501 if (_lower) {
1.502 for (EdgeIt e(_graph); e != INVALID; ++e)
1.503 - _flow[e] += (*_lower)[e];
1.504 + (*_flow)[e] += (*_lower)[e];
1.505 }
1.506 return true;
1.507 }